Bruno Crotman (crotman@gmail.com)
Trabalha há 6 anos na Empresa de Pesquisa Energética fomentando o trabalho com dados em todas as áreas da empresa
Trabalhou por 5 anos no Banco BBM/BBM Investimentos na área de Pesquisa Quantitativa
Doutorado (em curso): Engenharia de Software baseada em teoria dos jogos e ciência de dados, Unirio. Orientador: Prof. Márcio Barros
Mestrado: Engenharia de Software baseada em otimização, Unirio. Orientador: Prof. Márcio Barros
Graduação: Ciência da Computação, UFRJ
Introdução
Fundamentos da Linguagem
Introdução à Programação Funcional
Obtenção e Organização dos dados
Visualização de dados
Rudimentos de aprendizado estatístico
Comunicação
Séries Temporais ?
Práticas de Engenharia de Software ?
…da maldição do erro operacional
Fonte: (QLD 2016)
PS.: ele não é Nobel…
…e ainda…
Fonte: (Boddy 2016)
Meta final: que vários dos processos de análise de dados da empresa passem a ser feitos dentro do fluxo de trabalho do R.
Fonte:(Frigaard 2017)
Ao fim do curso o objetivo é que todos os fios da meada sejam puxados para que o aluno consiga continuar por si só usando a vasta documentação disponível.
Também amo o Excel, mas amo mais as seguintes vantagens:
Reprodutibilidade. Muito mais fácil refazer uma análise com código do que point and click
Menor risco operacional. A automatização é maior, a chance de erro na execução de um passo manual é nula
Menor risco de continuidade caso haja imprevistos com a equipe.
Maior flexibilidade. Virtualmente tudo é possível
Manutenção mais fácil
Controle de versão de forma profissional
Mais fácil do que parece
Ela é feita para lidar com dados
Comunidade de usuários gigante e cooperativa
Ferramentas poderosas para comunicação dos resultados, em documentos ou aplicações
Muitos pesquisadores em métodos quantitativos que estão no estado-da-arte publicam seus métodos em bibliotecas escritas em R
Prazeroso programar (Tidyverse, Shiny, R Markdown…)
Guerra de linguagens: R x Julia x Python. O mais importante é evitar Analysis Paralisis
Fonte:(Wickham and Grolemund 2016)
Fonte: (Frigaard 2017)
(Mello 2019)
É muito comum possuirmos dados gerados em planilhas ou em algum suporte de formato estruturado ou semi-estruturado.
Estes dados podem ser organizados de forma “tidy” para análise
Após a possível execução de modelos, podemos publicar os resultados.
Impacto da temperatura no consumo de energia elétrica
Localização de empresas e dutos
Estimador de posição de fundos multimercado
Cenários futuros para expansão da matriz energética
Localização dos equipamentos de geração distribuída
Equilíbrio de Nash pra vários cenários
Ajuda pra escolher o local de um consultório de psicologia
Formato de relatório. Análise mais direcionada e com conclusões do que exploratória
R é uma linguagem que é interpretada por um engine gratuito.
RStudio é o melhor ambiente de programação da linguagem R. A versão mais simples, que é totalmente funcional, é gratuita.
Na visualização padrão, ele oferece um console para execução de comandos e uma janela com a visualização dos environments, ou seja, das variáveis que ele guarda na sessão atual.
No console é possível executar comandos, como o que atribui valor a uma variável
Note que a atribuição é feita com <- e não com = como na maioria das linguagens.
Dica: o atalho alt + - gera o sinal de atribuição
Os comandos que não atribuem valor a uma variável são ecoados na tela
## [1] 3
Veja o [1] no console. O R considera que tudo é um vetor. É uma linguagem muito baseada em operações vetoriais. Isso facilita muito as coisas quando se lida com dados.
O console serve só para testes, aprendizado de novos comandos, debug, experiências etc.
Para executar adequadamente as atividades mais comuns de análise de dados, e para que elas sejam reprodutíveis, é necessária a criação de scripts.
Eles são salvos em um arquivo de extensão “.R”
Atalhos de teclado: ctrl+enter (rodar linhas selecionadas), ctrl+shift+enter (rodar script), ctrl+1 (foco no script), ctrl+2 (foco no console), ctrl+shift+F10 (reiniciar R), ctrl+shift+C (comentar/descomentar bloco) …
Refactoring
Document outline
Pane: Files/Plots/Packages/Help/Viewer
Pane: Environment/History/Connections/Git
Jobs
Controle de versão integrado com o Github
Cheat sheets
Verificação ortográfica (ainda não existe em portguês. Em inglês funciona bem)
Todo o material do curso está hospedado no Github, inclusive esta apresentação, escrita em RMarkdown.
Os exemplos de código, as imagens e os dados mostrados nesta apresentação estão inclusos no repositório do curso.
O repositório fica em github/crotman/cursoR.
Para baixar este repositório no RStudio, crie um projeto em File/New Project, do tipo Github e use o endereço do repositório: https://github.com/crotman/CursoR.git.
O projeto tanbém está disponível no RStudio Cloud: Projeto Curso R. A vantagem é que as bibliotecas já estão todas instaladas.
Todo material é disponibilizado sob a licença Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License
Para o R, simplificando para o escopo deste curso, as variáveis “armazenam” os seguintes tipos:
## [1] 1 2 3 4 5 6 7 8 9 10
## [[1]]
## [1] "oi"
##
## [[2]]
## [1] 1
## # A tibble: 10 x 2
## col1 col2
## <int> <int>
## 1 1 11
## 2 2 12
## 3 3 13
## 4 4 14
## 5 5 15
## 6 6 16
## 7 7 17
## 8 8 18
## 9 9 19
## 10 10 20
## [1] 3
## [1] "sou o b"
Existe orientação a objetos no R, mas não está no escopo deste curso
Note que não há variáveis que armazenam dado escalar, como já vimos.
Dentre os vetores há:
vetores atômicos (seus elementos são do mesmo tipo de dados primário)
listas (seus elementos, que são vetores atômicos, são de tipos de dados primários diferentes)
Tipos de vetores
Fonte: (Wickham 2019)
Tipos de dados primários:
Tipos primários
Fonte: (Wickham 2019)
## [1] FALSE
## [1] "integer"
## [1] "double"
## [1] 0.1
## [1] 1500
## [1] Inf
Uma das funções mais usadas do R é c(), que cria um vetor novo vetor combinando vetores.
## [1] 1 2 3
## [1] 1 2 3 4 5 6
## [1] 1.4 2.4 3.4 4.4 5.4 6.4 7.4 8.4 9.4
O operador : é usado para gerar um vetor com todos números que estão entre os operandos e são formados somando números inteiros ao primeiro operando.
## [1] 1 2 3 4 5 6 7 8 9 10
## [1] 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5
A função seq() é usada para criar um vetor de várias formas.
Numa das formas especifica-se o valor inicial, o valor final e o incremento entre elementos do vetor.
## [1] 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8
## [20] 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7
## [39] 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6
## [58] 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0 8.1 8.2 8.3 8.4 8.5
## [77] 8.6 8.7 8.8 8.9 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9
Note que chamamos a função sem especificar o nome dos parâmetros. Eles são recebidos pela função na ordem padrão definida pela função.
Mas no R também é possível passar parâmetros nomeados.
Clique em F1 enquanto tem o cursor em cima da função e veja a ordem dos parâmetros. Veja que outros parâmetros que não utilizamos. Podemos usar length.out ao invés de by:
## [1] 1 2 3 4 5 6 7 8 9 10
## [1] 1.00 3.25 5.50 7.75 10.00
Outro parâmetro, along.with, deixa que criemos um vetor num intervalo determinado e o mesmo número de elementos do vetor passado por este parâmetro.
## [1] 20.00000 28.88889 37.77778 46.66667 55.55556 64.44444 73.33333
## [8] 82.22222 91.11111 100.00000
Valores faltantes ou desconecidos são representados por NA
## [1] 1 NA
O valor NA quase sempre contamina os cálculos
## [1] NA
mas…
## [1] 1
A exceção são expressões que dão sempre o mesmo resultado independentemente do valor da variável
## [1] 1
## [1] TRUE
## [1] FALSE
A melhor forma de testar se existe um valor NA é is.na
## [1] FALSE TRUE FALSE
As operações do R são vetoriais. Numa operação entre um vetor e um escalar, a operação com o escalar é aplicada a cada elemento do vetor
## [1] 2 4 6 8 10
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Numa operação com vetores do mesmo tamanho, os elementos são pareados
## [1] 1 4 9 16 25 36 49 64 81 100
Um conceito importante é o de recycling.
Numa operação entre dois vetores de tamanhos diferentes, o vetor menor é repetido ciclicamente de forma a ficar com o mesmo tamanho do vetor maior.
Lembra que toda variável no R é um vetor?
Então… o escalar mostrado no primeiro código do slide anterior é um vetor de 1 elemento que sofre recycling
## [1] 1 4 3 8 5 12 7 16 9 20
Um vetor nomeado pode ser criado como abaixo
## primeiro segundo terceiro
## 1 2 3
Além de acessar o elemento do vetor pelo índice numérico (COMEÇANDO DE UM, pois R é uma linguagem 1-based ), é possível acessar os elementos pelos seus nomes
## primeiro
## 1
## primeiro
## 1
## terceiro
## 3
## terceiro
## 3
Os nomes dos elementos podem ser acessados com a função names()
## [1] "primeiro" "segundo" "terceiro"
Um vetor de nomes pode ser atribuído a names(). A mudança dos nomes é feita inline
## a b c
## 1 2 3
O acesso pode ser feito por um vetor, também.
## b c
## 2 3
As listas nomeadas funcionam de forma similar aos vetores nomeados.
## $primeiro
## [1] "a"
##
## $segundo
## [1] 2
##
## $terceiro
## [1] 1 2 3
Os elementos das listas são acessados com duplos colchetes
## [1] 1 2 3
## $um
## [1] "a"
##
## $dois
## [1] 2
##
## $tres
## [1] 1 2 3
Existem estruturas mais complexas construídas a partir de vetores e listas.
Data Frame
Matrix
Array
Factor
Estruturas que representam datas
Objetos (no paradigma de orientação a objetos)
Data Frame é a estrutura que mais vamos usar para nossas análises de dados.
Data Frame, e seu primo Tibble, são estruturas muito usadas em análises de dados feitas em R.
O Data Frame consiste em um conjunto de vetores/listas nomeados, com o mesmo número de elementos, que formam uma estrutura retangular, onde cada coluna é um vetor/lista e cada linha n contém o n-ésimo elemento dos vetores.
É similar, em muitas características, a uma tabela de banco de dados.
Essa estrutura é chave no paradigma “Tidy” que usaremos com as bibliotecas Tidyverse
Tibble é uma adaptação do Data Frame para análise de dados. Pra maioria dos efeitos eles têm o mesmo comportamento, mas os tibbles contêm algumas atualizações para os processos de análise de dados.
df <-
data.frame(
nome = c("João", "Maria", "Zezinho", "Juquinha"),
idade = c(7, 8, 9, 10),
altura = c(10, 11)
)
df## nome idade altura
## 1 João 7 10
## 2 Maria 8 11
## 3 Zezinho 9 10
## 4 Juquinha 10 11
#tibble não aceita recycling em vetores de tamanho diferente de 1
tib <-
#try evita que o erro paralise toda a execução do script
try(
tibble(
nome = c("João", "Maria", "Zezinho", "Juquinha"),
idade = c(7, 8, 9, 10),
altura = c(10, 11)
)
)## Error : Tibble columns must have compatible sizes.
## * Size 4: Existing data.
## * Size 2: Column `altura`.
## i Only values of size one are recycled.
Vamos usar bastante as list-columns. Ou seja, vendo o tibble como uma tabela, a coluna de um tibble pode ser um vetor (tipos primários iguais) ou uma lista (tipos primários diferentes).
Vamos usar list-columns por exemplo, ao aninhar conteúdos diferentes em cada linha de um tibble
Tibble é mais amigável a list-columns mas é possível usar list-columns em Data Frames também.
Tibble, por exemplo, aceita inicialização de uma list-column diretamente
## tibble [2 x 2] (S3: tbl_df/tbl/data.frame)
## $ a: num [1:2] 1 2
## $ b:List of 2
## ..$ : chr "a"
## ..$ : num 2
Data Frame se enrola…
## 'data.frame': 2 obs. of 3 variables:
## $ a : num 1 2
## $ b..a.: chr "a" "a"
## $ b.2 : num 2 2
Apesar de não poder ser inicializado desta forma, o Data Frame aceita list-column
## 'data.frame': 2 obs. of 2 variables:
## $ a: num 1 2
## $ b:List of 2
## ..$ : chr "a"
## ..$ : num 2
tibble com tribble()É possível criar um Tibble a partir de um código que parece uma tabela, ou seja, criar por linhas e não por colunas.
Isso é feito com a função tribble()
tribble(
~nome, ~idade, ~altura,
"Bruno", 41, 1.75,
"João", 23, 1.80,
"Maria", 29, 1.70,
"Zezinho", 31, 1.78
)## # A tibble: 4 x 3
## nome idade altura
## <chr> <dbl> <dbl>
## 1 Bruno 41 1.75
## 2 João 23 1.8
## 3 Maria 29 1.7
## 4 Zezinho 31 1.78
A linguagem oferece comandos de controle de fluxo similares aos de outras linguagens.
Podemos dividir os comandos de controle de fluxo em dois tipos:
choices: execução alternativa de comandos
loops: execução repetida de comandos
if, ifelseO comando if funciona para um valor lógico escalar
## [1] "2 mais 2 são 4"
Note o operador de comparação == e não =
A função if_else (da biblioteca dplyr) funciona de vetorial. if_else é mais rápida que a função ifelse da biblioteca base, mas só aceita argumentos de mesmo tipo no segundo e terceiro parâmetros
jogo_do_pim_silvio_santos <- if_else(
condition = 1:40 %% 4 == 0 ,
true = "PIM",
false = as.character(1:40)
)
jogo_do_pim_silvio_santos## [1] "1" "2" "3" "PIM" "5" "6" "7" "PIM" "9" "10" "11" "PIM"
## [13] "13" "14" "15" "PIM" "17" "18" "19" "PIM" "21" "22" "23" "PIM"
## [25] "25" "26" "27" "PIM" "29" "30" "31" "PIM" "33" "34" "35" "PIM"
## [37] "37" "38" "39" "PIM"
Note o operador %% e a função de coerção de tipo as.character
switch e case_whenA cláusula switch e a função dplyr::case_when evitam que o programador tenha que criar muitos if else aninhados
## [1] "começa com b"
Note que a condição vai sendo testada na ordem e stop gera um erro
case_when serve ao caso vetorial
## [1] "1" "par" "3" "par" "5" "par" "7" "par"
## [9] "9" "dezena" "11" "par" "13" "par" "15" "par"
## [17] "17" "par" "19" "dezena" "21" "par" "23" "par"
## [25] "25" "par" "27" "par" "29" "dezena" "31" "par"
## [33] "33" "par" "35" "par" "37" "par" "39" "dezena"
A cláusula de loop mais usada e mais versátil é for
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25
As cláusulas next e break modificam o comportamento, respectivamente caminhando direto para a próxima iteração e saindo do for
## [1] 1
## [1] 3
## [1] 5
## [1] 1
Vamos ver que quase sempre é desnecessário usar loop para as tarefas que vamos executar.
O caráter vetorial da linguagem, aliado a funcionalidades das bibliotecas, faz com que a grande maioria dos loops sejam desnecessários.
O código fica mais limpo e expressivo e mais rápido. Às vezes MUITO mais rápido. Isso ocorre por motivos além do escopo do curso (alocação de memória, código interpretado x código compilado em C++ etc.)
O código abaixo usa loop e programação funcional, respectivamente. Programação funcional será abordada posteriormente no material.
com_loop <- function(n){
x <- integer()
for (i in 1:n){
x <- c(x, i^2)
}
x
}
#programação funcional: aprenderemos posteriomente
sem_loop <- function(n){
x <- 1:n %>%
map_dbl(function(x){x^2})
x
}Abaixo as três formas de fazer a mesma conta que terão a performance avaliada
## [1] 1 4 9 16 25
## [1] 1 4 9 16 25
## [1] 1 4 9 16 25
A biblioteca bench oferece funções ótimas para avaliar a performance de pedaços pequenos de código.
resultados_perf <- mark(
sem_loop(1e4),
com_loop(1e4),
(1:1e4)^2
)
#aprenderemos o que é %>% e select() posteriormente
resultados_perf %>%
select(expression, min, median, `itr/sec` )## # A tibble: 3 x 4
## expression min median `itr/sec`
## <bch:expr> <bch:tm> <bch:tm> <dbl>
## 1 sem_loop(10000) 7.1ms 7.88ms 119.
## 2 com_loop(10000) 113.6ms 113.7ms 5.16
## 3 (1:10000)^2 14.6us 15.8us 16130.
Qual o valor de v1 ?
Qual o valor de v2
Complete a lacuna:
Qual o valor de x ?
Complete a lacuna para imprimir os quadrados dos números de 1 a 10
Monty Hall era uma espécie de Sílvio Santos juvenil (sub 80) americano.
Um dos seus jogos consistia em mostrar três portas ao otár… (ops) convidado. Atrás de uma delas sempre tem um carro. Atrás das outras, alguma brincadeira
Após a escolha inicial, o apresentador revelava uma das portas e perguntava se o convidado queria trocar a porta escolhida originalmente.
O que vocês acham? Melhor trocar, manter a escolha original ou tanto faz?
Note o que há de interessante no código (comentado):
A função sample()
A seleção negativa de elementos com -
A função replicate()
set.seed(88)
joga_monty_hall <- function(troca){
portas <- 1:3
#sample() sorteia elementos com ou sem reposição
porta_carro <- sample(portas, size = 1, replace = FALSE)
primeira_escolha <- 1
#Seleção negativa (retirando elementos)
portas_pra_revelar <- portas[-c(porta_carro, primeira_escolha)]
porta_revelada <- sample( c(portas_pra_revelar, portas_pra_revelar ), 1)
if(troca){
escolha <- portas[-c(primeira_escolha, porta_revelada)]
}
else{
escolha <- primeira_escolha
}
escolha == porta_carro
}
n <- 1000
#replicate executa múltiplas vezes um comando e armazena os resultados em uma estrutura única
troca <- replicate(n = n, joga_monty_hall(troca = TRUE))
fica <- replicate(n = n, joga_monty_hall(troca = FALSE))Resultados:
## [1] 0.675
## [1] 0.323
Vamos ver… mas dá pra simular sem saber quase nada.
Vamos usar uma das funções da família r<familia de distribuição de prob>(). Neste caso, a rbinom, que simula a distribução binomial, aquela que equivale a jogar n moedas para cima e ver quantas deram cara, ou seja, quantos desfechos chamados de “sucessos” são obtidos em n eventos com desfecho binário. Nestes eventos há uma probabilidade fixa para um desfecho considerado “sucesso.”
n_simul <- 10000
n_questoes <- 240
min_aprovacao <- 0.2
n_aprovado <- 240 * min_aprovacao
prob_questao <- 0.2
acertos <- rbinom(n = n_simul, size = n_questoes, prob = prob_questao )
sum(acertos >= n_aprovado)/n_simul ## [1] 0.5205
A chance é praticamente nula.
Na verdade, a grande massa da distribuição fica muito distante.
dado <- enframe(acertos/n_questoes)
mostra_chances <- function(acertos, n_questoes){
ggplot(enframe(acertos/n_questoes)) +
geom_density( aes(x = value)) +
scale_x_continuous(
labels = percent_format(accuracy = 1),
limits = c(0,1),
breaks = seq(0, 1, 0.1)
) +
labs(x ="% Acertos") +
geom_vline(xintercept = min_aprovacao, color = "red") +
theme_light()
}
mostra_chances(acertos, n_questoes)O exemplo anterior era muito simplista: ninguém chuta tudo.
Imagine que sabemos qual a chance de aparecer uma pergunta onde podemos descartar 0 alternativas, a chance de uma onde descartamos 1 e assim por diante.
Na função rbinom(), a cada simulação, sorteamos uma de dois tipos de desfecho para cada um dos n eventos. A rmultinom() sorteia um de k tipos de desfecho para cada um dos n eventos. Um vetor de k elementos indica a probabilidade de cada tipo de desfecho.
A rmultinom() nos ajuda a compor provas de acordo com a proficiência de quem vai fazer a prova.
#definindo a chance podermos eliminar 0, 1, 2, ... 4 alternativas
fracao_eliminar_questoes <- c( 0.1, 0.1, 0.2, 0.25 , 0.35 )
#definindo o número de questões
n_questoes_cada_elimina <- t(rmultinom(n_simul, size = n_questoes, fracao_eliminar_questoes))
probs_quando_elimina <- 1/(5:1)
acertos_concatenados <-
rbinom(
n = n_simul * 5 ,
size = as.vector(t(n_questoes_cada_elimina)),
prob = probs_quando_elimina
)## [,1] [,2] [,3] [,4] [,5]
## [1,] 23 24 50 47 96
## [2,] 24 17 48 64 87
## [3,] 23 30 43 56 88
## [4,] 20 25 51 52 92
## [1] 4 6 17 25 96 4 3 13 34 87 6 5 13 36 88 4 9 17 25 92
Como geramos os acertos concatenados em um vetor, o ideal é transformá-los numa matriz onde cada linha corresponde a uma simulação e cada coluna corresponde aos acertos conseguidos em cada tipo de questão. Usamos a função A matrix(), informando que estamos passando a matriz linha a linha e informando a quantidade de linhas.
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4 6 17 25 96
## [2,] 4 3 13 34 87
## [3,] 6 5 13 36 88
## [4,] 4 9 17 25 92
## [5,] 6 5 11 30 92
A rowSums() nos ajuda a ter a quantidade de acertos em cada simulação de prova.
## [1] 1
O fato de que funções são objetos de primeira classe no R, ou seja, objetos que têm propriedades iguais às de qualquer outro, possibilita a programação no estilo funcional.
A programação funcional inclui as formas de interação entre vetores e funções a partir de uma função
Neste curso vamos nos concentrar nos Functionals, funções que recebem função como argumento e devolvem um vetor ou uma lista
A função mais fundamental é map().
Ela recebe uma função e um vetor (ou lista) de \(n\) elementos
A função é chamada para cada elemento do vetor (ou lista), \(n\) vezes.
Os resultados da aplicação destas \(n\) execuções são devolvidos em uma lista de \(n\) elementos
Uma função é um objeto como outro qualquer e pode ser colocado em uma variável
map devolve uma lista com o resultado da execução de map em cada elemento
## [[1]]
## [1] 1
##
## [[2]]
## [1] 4
##
## [[3]]
## [1] 9
##
## [[4]]
## [1] 16
##
## [[5]]
## [1] 25
##
## [[6]]
## [1] 36
##
## [[7]]
## [1] 49
##
## [[8]]
## [1] 64
##
## [[9]]
## [1] 81
##
## [[10]]
## [1] 100
map pode retornar outros tipos de objeto## [1] 1 4 9 16 25 36 49 64 81 100
## [1] "1.000000" "4.000000" "9.000000" "16.000000" "25.000000"
## [6] "36.000000" "49.000000" "64.000000" "81.000000" "100.000000"
## # A tibble: 10 x 1
## x
## <dbl>
## 1 1
## 2 4
## 3 9
## 4 16
## 5 25
## 6 36
## 7 49
## 8 64
## 9 81
## 10 100
Funções podem ser declaradas inline.
## [[1]]
## [1] 1.0010565 1.0035960 0.9920089 0.9920322
##
## [[2]]
## [1] 2.006613 1.985193 1.995651 1.995989
##
## [[3]]
## [1] 2.986562 2.986561 2.987338 3.000187
##
## [[4]]
## [1] 3.977057 4.010954 4.005156 3.991333
##
## [[5]]
## [1] 5.008841 4.988631 4.997555 5.007538
Ou então shortcuts para funções podem ser usados
## [[1]]
## [1] 1.008859 1.012186 1.002763 1.004866
##
## [[2]]
## [1] 1.998966 1.992176 2.008423 2.011253
##
## [[3]]
## [1] 2.984938 3.009937 3.004963 3.008983
##
## [[4]]
## [1] 3.987746 3.985322 4.007886 3.991292
##
## [[5]]
## [1] 4.989955 4.996758 4.999391 5.015683
Os argumentos são repassados pela map para a função executada. Podemos então passar os argumentos para a map.
## [[1]]
## [1] 1.0081897 1.0112569 1.0143750 0.9860779
##
## [[2]]
## [1] 2.001435 1.977795 1.995766 2.000849
##
## [[3]]
## [1] 3.002678 2.983602 2.989823 2.994209
##
## [[4]]
## [1] 3.993385 3.994301 3.998010 3.996389
##
## [[5]]
## [1] 4.988824 4.996532 4.981557 4.991112
map() para dois argumentosA função map()possui versões para mais de um argumento.
map2 e suas modificações para retornos em outros tipos, como map2_dbl(), são preparadas para funções de dois argumentos.
tib_2_arg <- tribble(
~mean, ~sd,
0, 2,
0, 4,
1, 2,
1, 4,
)
map2(
.x = tib_2_arg$mean,
.y = tib_2_arg$sd,
.f = ~rnorm(n = 4, mean = .x, sd = .y)
)## [[1]]
## [1] 2.821674 -2.487257 -1.100323 -4.174959
##
## [[2]]
## [1] -1.1384564 -2.7152693 0.1741787 1.8220994
##
## [[3]]
## [1] -0.5682192 3.7054536 3.7375878 -1.3283584
##
## [[4]]
## [1] 3.417435 4.689878 7.298377 -1.680979
map() para mais de dois argumentospmap e suas modificações para retornos em outros tipos, como pmap_dbl(), são preparadas um número ilimitado de argumentos.
Uma lista de vetores nomeados deve ser passada para a pmap(). Esses são os argumentos passados para a função.
Lembre que um Data Frame/Tibble é exatamente isso: uma lista de vetores nomeados.
O detalhe é que os nomes dos vetores deve bater com o nome dos parâmetros da função
tib_n_arg <- tribble(
~mean, ~sd, ~n,
0, 2, 5,
0, 4, 5,
1, 2, 5,
1, 4, 5,
0, 2, 10,
0, 4, 10,
1, 2, 10,
1, 4, 10
)
pmap(.l = tib_n_arg, .f = rnorm)## [[1]]
## [1] 0.8674085 -4.5345372 -0.9266187 2.8163247 -1.0953331
##
## [[2]]
## [1] 5.136385 -3.057098 -2.295600 -1.794702 -1.311037
##
## [[3]]
## [1] 1.532084 -7.068620 2.369469 -1.790043 -2.132151
##
## [[4]]
## [1] 5.5805530 0.8476443 1.8301349 3.0455240 -6.1892686
##
## [[5]]
## [1] 2.5492313 0.7512555 -2.1446984 3.1894252 -3.5840746 2.2253177
## [7] -0.6031117 0.6509431 0.6480136 2.1763102
##
## [[6]]
## [1] 1.2978069 -1.7907063 -3.2601022 4.0173588 -2.3454343 -1.0336582
## [7] 0.6735530 -2.6442484 0.5021318 -2.2966865
##
## [[7]]
## [1] 3.2827845 0.8395007 -0.7688994 0.3239543 1.9593301 4.4368960
## [7] -3.6475308 1.1799538 1.9458926 3.2373521
##
## [[8]]
## [1] 2.330355 5.263333 1.419372 0.371877 -6.420001 5.127040 -5.708234
## [8] 2.881843 2.085133 3.274547
A programação funcional deixa o código mais expressivo e maioria das vezes mais performático
Existem algumas características avançadas da programação funcional que facilitam bastante o trabalho.
A principal delas é a facilidade com que podemos converter um código que roda em série para um código que roda em paralelo com a biblioteca furrr, que veremos adiante.
Dentro de todo o hype envolvendo Data Science, surgem as buzz words mais mirabolantes: machine learning, AI, Deep Learning…
Tudo isso é legal, mas a habilidade de preparar os dados para os modelos, preparar os hiperparâmetros e especificações alternativas ainda melhora muito a análise. Anote mais uma buzz word: FEATURE ENGINEERING
A agilidade de se tentar abordagens alterativas com os dados cresce muito quando dominamos a arte de manipular os data frames. Por isso o peso grande dado neste curso.
Arrumar os dados de forma que as linhas sejam eventos e as colunas sejam atributos do evento ajuda muito a rodar modelos e construir visualizações eficientemente. Esta forma de organização foi chamada de Tidy data por Hadley Wickham, o criador do conjunto de bibliotecas Tidyverse
O que é o evento e o que é o atributo pode variar até para diferentes usos do mesmo dado. Mas a prática ajuda a determinar isso.
%>%)O operador pipe, representado por %>% é extremamente útil para análise de dados no R com uso das bibliotecas Tidyverse
Normalmente os tratamentos de dados são feitos em múltiplos passos encadeados:
## # A tibble: 6 x 6
## country continent year lifeExp pop gdpPercap
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Afghanistan Asia 1952 28.8 8425333 779.
## 2 Afghanistan Asia 1957 30.3 9240934 821.
## 3 Afghanistan Asia 1962 32.0 10267083 853.
## 4 Afghanistan Asia 1967 34.0 11537966 836.
## 5 Afghanistan Asia 1972 36.1 13079460 740.
## 6 Afghanistan Asia 1977 38.4 14880372 786.
Vamos imaginar que queremos a média de PIB per capita por continente em 2007.
Note quanto código desnecessário há nestas linhas: variáveis que não precisavam ser nomeadas nem passadas explicitamente como parâmetro.
Este código desnecessário causa fadiga no programador, confunde o próprio autor do código e confunde ainda mais o leitor posterior do código.
#vamos cobrir essas funções de tratamento posteriormente
gapminder_07 <- filter(gapminder, year == 2007)
gapminder_07_group_continente <- group_by(gapminder_07, continent)
gapminder_media_gdp_continente <- summarise(
gapminder_07_group_continente, media_gdp = sum(gdpPercap * pop)/sum(pop)
)
resultado <- arrange(gapminder_media_gdp_continente, desc(media_gdp))
resultado## # A tibble: 5 x 2
## continent media_gdp
## <fct> <dbl>
## 1 Oceania 32885.
## 2 Europe 25244.
## 3 Americas 21603.
## 4 Asia 5432.
## 5 Africa 2561.
%>%) (cont.)O operador pipe %>% faz o seguinte:
x %>% y(z) = y(x,z)
Ou seja, o primeiro operando é enfiado como primeiro parâmetro da função que está no segundo operando.
Isso faz com que possamos escrever o código anterior assim:
resultado <- gapminder %>%
filter(year == 2007) %>%
group_by(continent) %>%
summarise(
media_gdp = sum(gdpPercap * pop) / sum(pop)
) %>%
arrange(desc(media_gdp))## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 5 x 2
## continent media_gdp
## <fct> <dbl>
## 1 Oceania 32885.
## 2 Europe 25244.
## 3 Americas 21603.
## 4 Asia 5432.
## 5 Africa 2561.
Note que agora podemos interpretar o código facilmente como uma série de comandos de tratamento em cima dos dados.
Não é por coincidência que as funções de tratamento das bibliotecas tidyverse que veremos adiante são verbos e recebem os dados como primeiro parâmetro.
Agora o mais importante de tudo: O ATALHO PARA O %>% É CTRL + SHIFT + M
dplyrdplyr é a biblioteca mais central do conjunto tidyverse
Fonte: (Courtiol 2019)
CRAN é o repositório de bibliotecas mantido pelo R com contribuição de populares.
Além de funcionalidades estatísticas e funcionalidades para lidar com dados, há dados e funcionalidades para buscar dados online.
Usaremos várias das bases como exemplo.
A primeira é a do Banco Mundial, muito rica para quem gosta de dados socioeconômicos: wbstats
Para acessar um indicador precisamos achá-lo na base de indicadores com a função wb_search()
#pattern é uma expressão regular. \\ serve para dizer que "(" é mesmo "("
#e não o ( usado nas operações de expressão regular (fora do escopo do curso)
indicadores <- wb_search(pattern = "GINI index \\(World Bank estimate\\)")
indicadores## # A tibble: 1 x 3
## indicator_id indicator indicator_desc
## <chr> <chr> <chr>
## 1 SI.POV.GINI GINI index (World B~ Gini index measures the extent to which the~
Sabendo o ID do indicador, podemos consultá-lo com a função wb_data()
#mrvnev é most recent non-empty values. Pode ser usado para buscar os n valores mais recentes
gini = wb_data(indicator = "SI.POV.GINI") %>%
filter(
!is.na(SI.POV.GINI)
)
head(gini)## # A tibble: 6 x 9
## iso2c iso3c country date SI.POV.GINI unit obs_status footnote last_updated
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <date>
## 1 AL ALB Albania 2017 33.2 <NA> <NA> Based on ~ 2020-12-16
## 2 AL ALB Albania 2016 33.7 <NA> <NA> Based on ~ 2020-12-16
## 3 AL ALB Albania 2015 32.9 <NA> <NA> Based on ~ 2020-12-16
## 4 AL ALB Albania 2014 34.6 <NA> <NA> Based on ~ 2020-12-16
## 5 AL ALB Albania 2012 29 <NA> <NA> Based on ~ 2020-12-16
## 6 AL ALB Albania 2008 30 <NA> <NA> Based on ~ 2020-12-16
A função select() é usada para selecionar colunas do data frame/tibble
## Rows: 1,649
## Columns: 9
## $ iso2c <chr> "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", ...
## $ iso3c <chr> "ALB", "ALB", "ALB", "ALB", "ALB", "ALB", "ALB", "ALB"...
## $ country <chr> "Albania", "Albania", "Albania", "Albania", "Albania",...
## $ date <dbl> 2017, 2016, 2015, 2014, 2012, 2008, 2005, 2002, 1996, ...
## $ SI.POV.GINI <dbl> 33.2, 33.7, 32.9, 34.6, 29.0, 30.0, 30.6, 31.7, 27.0, ...
## $ unit <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ obs_status <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ footnote <chr> "Based on data from HBS. Estimated from unit-record co...
## $ last_updated <date> 2020-12-16, 2020-12-16, 2020-12-16, 2020-12-16, 2020-...
## # A tibble: 6 x 4
## country date SI.POV.GINI iso3c
## <chr> <dbl> <dbl> <chr>
## 1 Albania 2017 33.2 ALB
## 2 Albania 2016 33.7 ALB
## 3 Albania 2015 32.9 ALB
## 4 Albania 2014 34.6 ALB
## 5 Albania 2012 29 ALB
## 6 Albania 2008 30 ALB
É possível usar a seleção negativa assim como fizemos com vetores
## # A tibble: 6 x 3
## country date SI.POV.GINI
## <chr> <dbl> <dbl>
## 1 Albania 2017 33.2
## 2 Albania 2016 33.7
## 3 Albania 2015 32.9
## 4 Albania 2014 34.6
## 5 Albania 2012 29
## 6 Albania 2008 30
Algumas funções helpers da dplyr nos ajudam a usar a função select e são muito úteis para tratamentos mais elaborados.
Pra mostrar mais funcionalidades da função select, vamos usar uma base com dados eleitorais brasileiros, que retorna mais colunas
## Rows: 29,146
## Columns: 58
## $ DATA_GERACAO <chr> "01/10/2020", "01/10/2020", "01/10/2...
## $ HORA_GERACAO <time> 11:48:48, 11:48:48, 11:48:48, 11:48...
## $ ANO_ELEICAO <dbl> 2018, 2018, 2018, 2018, 2018, 2018, ...
## $ COD_TIPO_ELEICAO <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...
## $ NOME_TIPO_ELEICAO <chr> "ELEIÇÃO ORDINÁRIA", "ELEIÇÃO ORDINÁ...
## $ NUM_TURNO <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ COD_ELEICAO <dbl> 297, 297, 297, 297, 297, 297, 297, 2...
## $ DESCRICAO_ELEICAO <chr> "Eleições Gerais Estaduais 2018", "E...
## $ DATA_ELEICAO <chr> "07/10/2018", "07/10/2018", "07/10/2...
## $ ABRANGENCIA <chr> "ESTADUAL", "ESTADUAL", "ESTADUAL", ...
## $ SIGLA_UF <chr> "AC", "AC", "AC", "AC", "AC", "AC", ...
## $ SIGLA_UE <chr> "AC", "AC", "AC", "AC", "AC", "AC", ...
## $ DESCRICAO_UE <chr> "ACRE", "ACRE", "ACRE", "ACRE", "ACR...
## $ CODIGO_CARGO <dbl> 7, 6, 7, 6, 6, 7, 7, 7, 7, 7, 7, 7, ...
## $ DESCRICAO_CARGO <chr> "DEPUTADO ESTADUAL", "DEPUTADO FEDER...
## $ SEQUENCIAL_CANDIDATO <dbl> 10000610652, 10000608406, 1000060230...
## $ NUMERO_CANDIDATO <dbl> 11111, 3131, 12300, 5555, 2722, 1778...
## $ NOME_CANDIDATO <chr> "GEHLEN DINIZ ANDRADE", "JEFFERSON B...
## $ NOME_URNA_CANDIDATO <chr> "GERLEN DINIZ", "JEFFERSON BARROSO "...
## $ NOME_SOCIAL_CANDIDATO <chr> "#NULO#", "#NULO#", "#NULO#", "#NULO...
## $ CPF_CANDIDATO <chr> "35954590249", "85546976268", "30843...
## $ EMAIL_CANDIDATO <chr> "ANGELMFERREIRA@GMAIL.COM", "JEFFERS...
## $ COD_SITUACAO_CANDIDATURA <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, ...
## $ DES_SITUACAO_CANDIDATURA <chr> "APTO", "APTO", "APTO", "APTO", "APT...
## $ COD_DETALHE_SITUACAO_CAND <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ...
## $ DES_DETALHE_SITUACAO_CAND <chr> "DEFERIDO", "DEFERIDO", "DEFERIDO", ...
## $ TIPO_AGREMIACAO <chr> "COLIGAÇÃO", "COLIGAÇÃO", "PARTIDO I...
## $ NUMERO_PARTIDO <dbl> 11, 31, 12, 55, 27, 17, 17, 19, 23, ...
## $ SIGLA_PARTIDO <chr> "PP", "PHS", "PDT", "PSD", "DC", "PS...
## $ NOME_PARTIDO <chr> "PROGRESSISTAS", "PARTIDO HUMANISTA ...
## $ CODIGO_LEGENDA <dbl> 10000050286, 10000050220, 1000005006...
## $ NOME_COLIGACAO <chr> "MUDANÇA COM COMPROMISSO", "Frente P...
## $ COMPOSICAO_LEGENDA <chr> "PP / PMN / PPS / PTC / PR", "PT / P...
## $ CODIGO_NACIONALIDADE <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ DESCRICAO_NACIONALIDADE <chr> "BRASILEIRA NATA", "BRASILEIRA NATA"...
## $ SIGLA_UF_NASCIMENTO <chr> "AC", "AC", "AC", "AC", "AC", "AC", ...
## $ CODIGO_MUNICIPIO_NASCIMENTO <dbl> -3, -3, -3, -3, -3, -3, -3, -3, -3, ...
## $ NOME_MUNICIPIO_NASCIMENTO <chr> "RIO BRANCO", "RIO BRANCO ", "FEIJÓ"...
## $ DATA_NASCIMENTO <chr> "08/04/1974", "11/11/1994", "23/03/1...
## $ IDADE_DATA_POSSE <dbl> 46, 26, 52, 56, 41, 61, 37, 48, 44, ...
## $ NUM_TITULO_ELEITORAL_CANDIDATO <chr> "002092722470", "006267492410", "000...
## $ CODIGO_SEXO <dbl> 2, 2, 2, 2, 4, 2, 2, 4, 2, 4, 2, 2, ...
## $ DESCRICAO_SEXO <chr> "MASCULINO", "MASCULINO", "MASCULINO...
## $ COD_GRAU_INSTRUCAO <dbl> 8, 8, 6, 8, 6, 2, 8, 4, 6, 6, 6, 3, ...
## $ DESCRICAO_GRAU_INSTRUCAO <chr> "SUPERIOR COMPLETO", "SUPERIOR COMPL...
## $ CODIGO_ESTADO_CIVIL <dbl> 3, 1, 1, 9, 3, 3, 3, 3, 1, 3, 3, 1, ...
## $ DESCRICAO_ESTADO_CIVIL <chr> "CASADO(A)", "SOLTEIRO(A)", "SOLTEIR...
## $ CODIGO_COR_RACA <chr> "01", "01", "03", "01", "03", "03", ...
## $ DESCRICAO_COR_RACA <chr> "BRANCA", "BRANCA", "PARDA", "BRANCA...
## $ CODIGO_OCUPACAO <dbl> 277, 999, 257, 257, 257, 999, 999, 9...
## $ DESCRICAO_OCUPACAO <chr> "DEPUTADO", "OUTROS", "EMPRESÁRIO", ...
## $ DESPESA_MAX_CAMPANHA <dbl> 1000000, 2500000, 1000000, 2500000, ...
## $ COD_SIT_TOT_TURNO <dbl> 2, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, ...
## $ DESC_SIT_TOT_TURNO <chr> "ELEITO POR QP", "SUPLENTE", "SUPLEN...
## $ SITUACAO_REELEICAO <chr> "S", "N", "N", "N", "N", "N", "N", "...
## $ SITUACAO_DECLARAR_BENS <chr> "S", "N", "S", "S", "N", "N", "N", "...
## $ NUMERO_PROTOCOLO_CANDIDATURA <chr> "-1", "-1", "-1", "-1", "-1", "-1", ...
## $ NUMERO_PROCESSO <chr> "06006019120186010000", "06005594220...
candidatos_select <- candidatos %>%
select(ends_with("candidato"))
datatable(head(candidatos_select))A função helper num_range ajuda a encontrar colunas do tipo prefixo_n. Isso é muito comum em bases de dados
A biblioteca worldmet retorna dados de estações meteorológicas espalhadas pelo planeta
Primeiro é necessário encontrar o código da base desejada
A função abaixo retorna os dados de uma estação. Veja que alguns campos têm um sufixo _n
## Rows: 8,760
## Columns: 25
## $ code <fct> 037720-99999, 037720-99999, 037720-99999, 037720-99999,...
## $ station <fct> "HEATHROW, UK", "HEATHROW, UK", "HEATHROW, UK", "HEATHR...
## $ date <dttm> 2014-01-01 00:00:00, 2014-01-01 01:00:00, 2014-01-01 0...
## $ latitude <dbl> 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 5...
## $ longitude <dbl> -0.461389, -0.461389, -0.461389, -0.461389, -0.461389, ...
## $ elev <dbl> 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29,...
## $ ws <dbl> 6.000000, 5.866667, 4.433333, 3.266667, 3.600000, 4.433...
## $ wd <dbl> 200.3347, 190.0000, 180.0000, 166.8408, 172.9098, 180.0...
## $ air_temp <dbl> 7.233333, 6.600000, 6.100000, 6.766667, 6.833333, 7.333...
## $ atmos_pres <dbl> 1002.1, 1002.0, 1001.9, 1001.7, 1001.1, 1000.7, 999.9, ...
## $ visibility <dbl> 15000.000, 30000.000, 40000.000, 40000.000, 24000.000, ...
## $ dew_point <dbl> 4.833333, 4.433333, 4.866667, 4.966667, 5.100000, 5.733...
## $ RH <dbl> 84.93283, 86.22546, 91.92697, 88.45555, 88.87875, 89.74...
## $ ceil_hgt <dbl> 1030.0000, NA, NA, 732.0000, 975.0000, NA, 914.0000, 13...
## $ cl_1 <dbl> 1.666667, 1.666667, 2.666667, 3.333333, 1.666667, 3.000...
## $ cl_2 <dbl> 4.500000, NA, NA, 7.000000, 3.500000, NA, 4.000000, 6.6...
## $ cl_3 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ cl <dbl> 3.666667, 1.666667, 2.666667, 4.333333, 3.000000, 3.000...
## $ cl_1_height <dbl> 687.6667, 908.0000, 1150.0000, 1098.0000, 920.3333, 625...
## $ cl_2_height <dbl> 1030.0000, NA, NA, 732.0000, 1987.5000, NA, 914.0000, 1...
## $ cl_3_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ precip_12 <dbl> NA, NA, NA, NA, NA, NA, 8, NA, NA, NA, NA, NA, NA, NA, ...
## $ precip_6 <dbl> 6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20, NA, ...
## $ precip <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3.333333, 3.333...
## $ pwc <chr> " rain showers or intermittent rain, slight", " rain", ...
A função helper num_range ajuda a selecionar essas colunas com prefixo comum e um sufixo numérico
dados_heathrow_select <- dados_heathrow %>%
select(
date,
num_range("cl_", 1:3 ),
num_range("precip_", c(6, 12))
)
head(dados_heathrow_select)## # A tibble: 6 x 6
## date cl_1 cl_2 cl_3 precip_6 precip_12
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2014-01-01 00:00:00 1.67 4.5 NA 6 NA
## 2 2014-01-01 01:00:00 1.67 NA NA NA NA
## 3 2014-01-01 02:00:00 2.67 NA NA NA NA
## 4 2014-01-01 03:00:00 3.33 7 NA NA NA
## 5 2014-01-01 04:00:00 1.67 3.5 NA NA NA
## 6 2014-01-01 05:00:00 3 NA NA NA NA
Outra função útil é a everything(), que ajuda, por exemplo, a passar algumas colunas para o início do tibble.
dados_heathrow_select <- dados_heathrow %>%
select(
date,
air_temp,
everything()
)
glimpse(dados_heathrow_select)## Rows: 8,760
## Columns: 25
## $ date <dttm> 2014-01-01 00:00:00, 2014-01-01 01:00:00, 2014-01-01 0...
## $ air_temp <dbl> 7.233333, 6.600000, 6.100000, 6.766667, 6.833333, 7.333...
## $ code <fct> 037720-99999, 037720-99999, 037720-99999, 037720-99999,...
## $ station <fct> "HEATHROW, UK", "HEATHROW, UK", "HEATHROW, UK", "HEATHR...
## $ latitude <dbl> 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 5...
## $ longitude <dbl> -0.461389, -0.461389, -0.461389, -0.461389, -0.461389, ...
## $ elev <dbl> 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29,...
## $ ws <dbl> 6.000000, 5.866667, 4.433333, 3.266667, 3.600000, 4.433...
## $ wd <dbl> 200.3347, 190.0000, 180.0000, 166.8408, 172.9098, 180.0...
## $ atmos_pres <dbl> 1002.1, 1002.0, 1001.9, 1001.7, 1001.1, 1000.7, 999.9, ...
## $ visibility <dbl> 15000.000, 30000.000, 40000.000, 40000.000, 24000.000, ...
## $ dew_point <dbl> 4.833333, 4.433333, 4.866667, 4.966667, 5.100000, 5.733...
## $ RH <dbl> 84.93283, 86.22546, 91.92697, 88.45555, 88.87875, 89.74...
## $ ceil_hgt <dbl> 1030.0000, NA, NA, 732.0000, 975.0000, NA, 914.0000, 13...
## $ cl_1 <dbl> 1.666667, 1.666667, 2.666667, 3.333333, 1.666667, 3.000...
## $ cl_2 <dbl> 4.500000, NA, NA, 7.000000, 3.500000, NA, 4.000000, 6.6...
## $ cl_3 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ cl <dbl> 3.666667, 1.666667, 2.666667, 4.333333, 3.000000, 3.000...
## $ cl_1_height <dbl> 687.6667, 908.0000, 1150.0000, 1098.0000, 920.3333, 625...
## $ cl_2_height <dbl> 1030.0000, NA, NA, 732.0000, 1987.5000, NA, 914.0000, 1...
## $ cl_3_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ precip_12 <dbl> NA, NA, NA, NA, NA, NA, 8, NA, NA, NA, NA, NA, NA, NA, ...
## $ precip_6 <dbl> 6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20, NA, ...
## $ precip <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3.333333, 3.333...
## $ pwc <chr> " rain showers or intermittent rain, slight", " rain", ...
Outra forma mais nova de fazer isso é usando a função relocate(), que facilita a ordenação das colunas de um data frame
dados_heathrow_relocate <- dados_heathrow %>%
relocate(
date,
air_temp
)
glimpse(dados_heathrow_relocate)## Rows: 8,760
## Columns: 25
## $ date <dttm> 2014-01-01 00:00:00, 2014-01-01 01:00:00, 2014-01-01 0...
## $ air_temp <dbl> 7.233333, 6.600000, 6.100000, 6.766667, 6.833333, 7.333...
## $ code <fct> 037720-99999, 037720-99999, 037720-99999, 037720-99999,...
## $ station <fct> "HEATHROW, UK", "HEATHROW, UK", "HEATHROW, UK", "HEATHR...
## $ latitude <dbl> 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 5...
## $ longitude <dbl> -0.461389, -0.461389, -0.461389, -0.461389, -0.461389, ...
## $ elev <dbl> 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29,...
## $ ws <dbl> 6.000000, 5.866667, 4.433333, 3.266667, 3.600000, 4.433...
## $ wd <dbl> 200.3347, 190.0000, 180.0000, 166.8408, 172.9098, 180.0...
## $ atmos_pres <dbl> 1002.1, 1002.0, 1001.9, 1001.7, 1001.1, 1000.7, 999.9, ...
## $ visibility <dbl> 15000.000, 30000.000, 40000.000, 40000.000, 24000.000, ...
## $ dew_point <dbl> 4.833333, 4.433333, 4.866667, 4.966667, 5.100000, 5.733...
## $ RH <dbl> 84.93283, 86.22546, 91.92697, 88.45555, 88.87875, 89.74...
## $ ceil_hgt <dbl> 1030.0000, NA, NA, 732.0000, 975.0000, NA, 914.0000, 13...
## $ cl_1 <dbl> 1.666667, 1.666667, 2.666667, 3.333333, 1.666667, 3.000...
## $ cl_2 <dbl> 4.500000, NA, NA, 7.000000, 3.500000, NA, 4.000000, 6.6...
## $ cl_3 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ cl <dbl> 3.666667, 1.666667, 2.666667, 4.333333, 3.000000, 3.000...
## $ cl_1_height <dbl> 687.6667, 908.0000, 1150.0000, 1098.0000, 920.3333, 625...
## $ cl_2_height <dbl> 1030.0000, NA, NA, 732.0000, 1987.5000, NA, 914.0000, 1...
## $ cl_3_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ precip_12 <dbl> NA, NA, NA, NA, NA, NA, 8, NA, NA, NA, NA, NA, NA, NA, ...
## $ precip_6 <dbl> 6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20, NA, ...
## $ precip <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3.333333, 3.333...
## $ pwc <chr> " rain showers or intermittent rain, slight", " rain", ...
Por padrão, ela traz as colunas para a frente, mas é possível movê-las para antes ou depois de uma outra coluna, usando o parâmetro .before OU EXCLUSIVO .after
dados_heathrow_relocate <- dados_heathrow %>%
relocate(
date,
air_temp,
.before = station
)
glimpse(dados_heathrow_relocate)## Rows: 8,760
## Columns: 25
## $ code <fct> 037720-99999, 037720-99999, 037720-99999, 037720-99999,...
## $ date <dttm> 2014-01-01 00:00:00, 2014-01-01 01:00:00, 2014-01-01 0...
## $ air_temp <dbl> 7.233333, 6.600000, 6.100000, 6.766667, 6.833333, 7.333...
## $ station <fct> "HEATHROW, UK", "HEATHROW, UK", "HEATHROW, UK", "HEATHR...
## $ latitude <dbl> 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 51.4775, 5...
## $ longitude <dbl> -0.461389, -0.461389, -0.461389, -0.461389, -0.461389, ...
## $ elev <dbl> 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29, 25.29,...
## $ ws <dbl> 6.000000, 5.866667, 4.433333, 3.266667, 3.600000, 4.433...
## $ wd <dbl> 200.3347, 190.0000, 180.0000, 166.8408, 172.9098, 180.0...
## $ atmos_pres <dbl> 1002.1, 1002.0, 1001.9, 1001.7, 1001.1, 1000.7, 999.9, ...
## $ visibility <dbl> 15000.000, 30000.000, 40000.000, 40000.000, 24000.000, ...
## $ dew_point <dbl> 4.833333, 4.433333, 4.866667, 4.966667, 5.100000, 5.733...
## $ RH <dbl> 84.93283, 86.22546, 91.92697, 88.45555, 88.87875, 89.74...
## $ ceil_hgt <dbl> 1030.0000, NA, NA, 732.0000, 975.0000, NA, 914.0000, 13...
## $ cl_1 <dbl> 1.666667, 1.666667, 2.666667, 3.333333, 1.666667, 3.000...
## $ cl_2 <dbl> 4.500000, NA, NA, 7.000000, 3.500000, NA, 4.000000, 6.6...
## $ cl_3 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ cl <dbl> 3.666667, 1.666667, 2.666667, 4.333333, 3.000000, 3.000...
## $ cl_1_height <dbl> 687.6667, 908.0000, 1150.0000, 1098.0000, 920.3333, 625...
## $ cl_2_height <dbl> 1030.0000, NA, NA, 732.0000, 1987.5000, NA, 914.0000, 1...
## $ cl_3_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ precip_12 <dbl> NA, NA, NA, NA, NA, NA, 8, NA, NA, NA, NA, NA, NA, NA, ...
## $ precip_6 <dbl> 6, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20, NA, ...
## $ precip <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3.333333, 3.333...
## $ pwc <chr> " rain showers or intermittent rain, slight", " rain", ...
rename()rename() é usada para modificar os nomes das colunas. Ela renomeia as colunas indicadas e mantném as outras.
rename()## # A tibble: 6 x 25
## code estacao data latitude longitude elev ws wd
## <fct> <fct> <dttm> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0377~ HEATHR~ 2014-01-01 00:00:00 51.5 -0.461 25.3 6 200.
## 2 0377~ HEATHR~ 2014-01-01 01:00:00 51.5 -0.461 25.3 5.87 190
## 3 0377~ HEATHR~ 2014-01-01 02:00:00 51.5 -0.461 25.3 4.43 180
## 4 0377~ HEATHR~ 2014-01-01 03:00:00 51.5 -0.461 25.3 3.27 167.
## 5 0377~ HEATHR~ 2014-01-01 04:00:00 51.5 -0.461 25.3 3.6 173.
## 6 0377~ HEATHR~ 2014-01-01 05:00:00 51.5 -0.461 25.3 4.43 180
## # ... with 17 more variables: temperatura <dbl>, atmos_pres <dbl>,
## # visibility <dbl>, dew_point <dbl>, RH <dbl>, ceil_hgt <dbl>, cl_1 <dbl>,
## # cl_2 <dbl>, cl_3 <dbl>, cl <dbl>, cl_1_height <dbl>, cl_2_height <dbl>,
## # cl_3_height <dbl>, precip_12 <dbl>, precip_6 <dbl>, precip <dbl>, pwc <chr>
rename() x select()select() também pode ser usado para renomear colunas, mas mantém apenas as colunas citadas
## # A tibble: 6 x 3
## data estacao temperatura
## <dttm> <fct> <dbl>
## 1 2014-01-01 00:00:00 HEATHROW, UK 7.23
## 2 2014-01-01 01:00:00 HEATHROW, UK 6.6
## 3 2014-01-01 02:00:00 HEATHROW, UK 6.1
## 4 2014-01-01 03:00:00 HEATHROW, UK 6.77
## 5 2014-01-01 04:00:00 HEATHROW, UK 6.83
## 6 2014-01-01 05:00:00 HEATHROW, UK 7.33
rename_with()A função rename_with() possibilita que renomeemos várias colunas ao mesmo tempo com o uso de uma função.
Aqui fazemos uso de uma função da biblioteca stringr, que trata strings, ou seja, cadeias de caracteres.
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
Primeiro podemos nos incomodar com as letras maiúsculas e o “.” separando as palavras em vez de “_”
iris_padrao <- iris %>%
rename_with(
.fn = str_to_lower
) %>%
rename_with(
.fn = function(x){
str_replace(
string = x,
pattern = "\\.",
replacement = "_"
)}
)
str(iris_padrao)## 'data.frame': 150 obs. of 5 variables:
## $ sepal_length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ sepal_width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ petal_length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ petal_width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
Podemos também querer colocar um sufixo nos nomes dos campos com a unidade
Vamos usar aqui o “~” criando um atalho para a função, que vai receber cada nome de coluna como .x
A função str_glue() possibilita a criação de novas strings a partir de variáveis já existentes de uma forma expressiva.
iris_unidade <- iris_padrao %>%
rename_with(
.cols = -species ,
.fn = ~str_glue("{.x}_inches")
)
str(iris_unidade)## 'data.frame': 150 obs. of 5 variables:
## $ sepal_length_inches: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ sepal_width_inches : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ petal_length_inches: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ petal_width_inches : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
A função mutate() é usada para criar novas colunas no tibble, mantendo as outras.
mutate()Notando que a coluna DATA_ELEICAO é um caracter, vamos criar uma coluna de tipo data.
## [1] "character"
O jeito mais fácil de fazer isso é usando uma das funções da biblioteca lubridate
candidatos_com_data <- candidatos %>%
select(1:10) %>%
mutate(DATA_ELEICAO_TIPO_DATA = dmy(DATA_ELEICAO)) %>%
select(DATA_ELEICAO, DATA_ELEICAO_TIPO_DATA)
head(candidatos_com_data)## # A tibble: 6 x 2
## DATA_ELEICAO DATA_ELEICAO_TIPO_DATA
## <chr> <date>
## 1 07/10/2018 2018-10-07
## 2 07/10/2018 2018-10-07
## 3 07/10/2018 2018-10-07
## 4 07/10/2018 2018-10-07
## 5 07/10/2018 2018-10-07
## 6 07/10/2018 2018-10-07
É possível substituir um campo existente
candidatos_com_data <- candidatos %>%
select(1:10) %>%
mutate(DATA_ELEICAO = dmy(DATA_ELEICAO)) %>%
select(DATA_ELEICAO)
head(candidatos_com_data)## # A tibble: 6 x 1
## DATA_ELEICAO
## <date>
## 1 2018-10-07
## 2 2018-10-07
## 3 2018-10-07
## 4 2018-10-07
## 5 2018-10-07
## 6 2018-10-07
mutate()funções derivadas da mutate possibilitam a alteração de várias colunas ao mesmo tempo, usando os mesmos helpers que já vimos para a select e uma função à escolha.
Veremos uma forma mais nova de fazer isso, com a função across(), mas essa forma continua a funcionar e é usada em muitos códigos existentes.
candidatos_com_data <- candidatos %>%
select(1:10) %>%
mutate_at(vars(starts_with("DATA_")), dmy ) %>%
select(starts_with("DATA_"))
head(candidatos_com_data)## # A tibble: 6 x 2
## DATA_GERACAO DATA_ELEICAO
## <date> <date>
## 1 2020-10-01 2018-10-07
## 2 2020-10-01 2018-10-07
## 3 2020-10-01 2018-10-07
## 4 2020-10-01 2018-10-07
## 5 2020-10-01 2018-10-07
## 6 2020-10-01 2018-10-07
mutate() (cont.):Outras funções úteis são as que fazem operações acumuladas e as operações de lag() e lead().
BETS é uma biblioteca criada pela FGV que dá acesso a séries temporais econômicas
## # A tibble: 33 x 7
## code description unit periodicity start last_value source
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1 Exchange rate - Free - U~ c.m.u~ D 28/11~ 26/03/2018 Sisbace~
## 2 10813 Exchange rate - Free - U~ c.m.u~ D 28/11~ 03/05/2018 Sisbace~
## 3 11753 Real effective exchange ~ Index M 31/01~ mar/2018 BCB-DST~
## 4 11758 Real effective exchange ~ Index M 31/01~ mar/2018 BCB-DST~
## 5 11763 Real effective exchange ~ Index M 31/01~ mar/2018 BCB-DST~
## 6 11768 Real effective exchange ~ Index M 31/01~ mar/2018 BCB-DST~
## 7 17887 Exchange rate (period av~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## 8 17888 Exchange rate (period av~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## 9 18441 Exchange rate (end of pe~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## 10 18442 Exchange rate (end of pe~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## # ... with 23 more rows
No código abaixo, calculamos o retorno da série, a volatilidade histórica e a volatilidade EWMA
dolar <- BETS::BETSget(1)
dolar_com_info <- dolar %>%
filter(date > ymd("1994-07-01")) %>%
arrange(date) %>%
mutate(
retorno = (value - lag(value))/value,
retorno_quad = retorno^2,
dia = row_number(),
fator_ewma = (1/0.94)^dia*1e-20
) %>%
filter(!is.na(retorno))
dolar_com_vol <- dolar_com_info %>%
mutate(vol = sqrt(cumvar(retorno)) * sqrt(252) ) %>%
mutate(vol_ewma = sqrt(cumsum(retorno_quad * fator_ewma)/cumsum(fator_ewma)) * sqrt(252) ) %>%
rename(dolar = value) %>%
select(
date,
dolar,
retorno,
vol,
vol_ewma
)
datatable(dolar_com_vol) %>%
formatPercentage(c("retorno", "vol", "vol_ewma"), 2)dolar_ajeitado <- dolar_com_vol %>%
pivot_longer(cols = -date, names_to = "variavel", values_to = "valor")
dolar_ajeitado %>%
ggplot() +
geom_line(aes(x = date, y = valor)) +
facet_grid( variavel ~ . , scales = "free") +
theme_light() transmute():transmute() cria colunas e mantém apenas as colunas citadas
## # A tibble: 6 x 3
## ano pais pib
## <int> <fct> <dbl>
## 1 1952 Afghanistan 6567086330.
## 2 1957 Afghanistan 7585448670.
## 3 1962 Afghanistan 8758855797.
## 4 1967 Afghanistan 9648014150.
## 5 1972 Afghanistan 9678553274.
## 6 1977 Afghanistan 11697659231.
A função arrange serve para ordenar as linhas do tibble.
arrange():## # A tibble: 6 x 25
## code station date latitude longitude elev ws wd
## <fct> <fct> <dttm> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0377~ HEATHR~ 2014-01-01 00:00:00 51.5 -0.461 25.3 6 200.
## 2 0377~ HEATHR~ 2014-01-01 01:00:00 51.5 -0.461 25.3 5.87 190
## 3 0377~ HEATHR~ 2014-01-01 02:00:00 51.5 -0.461 25.3 4.43 180
## 4 0377~ HEATHR~ 2014-01-01 03:00:00 51.5 -0.461 25.3 3.27 167.
## 5 0377~ HEATHR~ 2014-01-01 04:00:00 51.5 -0.461 25.3 3.6 173.
## 6 0377~ HEATHR~ 2014-01-01 05:00:00 51.5 -0.461 25.3 4.43 180
## # ... with 17 more variables: air_temp <dbl>, atmos_pres <dbl>,
## # visibility <dbl>, dew_point <dbl>, RH <dbl>, ceil_hgt <dbl>, cl_1 <dbl>,
## # cl_2 <dbl>, cl_3 <dbl>, cl <dbl>, cl_1_height <dbl>, cl_2_height <dbl>,
## # cl_3_height <dbl>, precip_12 <dbl>, precip_6 <dbl>, precip <dbl>, pwc <chr>
A função desc() permite a ordenação decrescente
## # A tibble: 6 x 25
## code station date latitude longitude elev ws wd
## <fct> <fct> <dttm> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0377~ HEATHR~ 2014-12-31 23:00:00 51.5 -0.461 25.3 4.27 204.
## 2 0377~ HEATHR~ 2014-12-31 22:00:00 51.5 -0.461 25.3 4.77 200
## 3 0377~ HEATHR~ 2014-12-31 21:00:00 51.5 -0.461 25.3 5.5 197.
## 4 0377~ HEATHR~ 2014-12-31 20:00:00 51.5 -0.461 25.3 4.8 197.
## 5 0377~ HEATHR~ 2014-12-31 19:00:00 51.5 -0.461 25.3 4.43 183.
## 6 0377~ HEATHR~ 2014-12-31 18:00:00 51.5 -0.461 25.3 4.6 193.
## # ... with 17 more variables: air_temp <dbl>, atmos_pres <dbl>,
## # visibility <dbl>, dew_point <dbl>, RH <dbl>, ceil_hgt <dbl>, cl_1 <dbl>,
## # cl_2 <dbl>, cl_3 <dbl>, cl <dbl>, cl_1_height <dbl>, cl_2_height <dbl>,
## # cl_3_height <dbl>, precip_12 <dbl>, precip_6 <dbl>, precip <dbl>, pwc <chr>
filter():filter() seleciona colunas de acordo com os seus valores
filter() (cont.):Filtrando países (note o operador %in% )
gapminder %>%
filter(country %in% c("Brazil", "Argentina", "Chile")) %>%
ggplot() +
geom_line(aes(x = year, y = gdpPercap, color = country )) +
geom_point(aes(x = year, y = gdpPercap, color = country )) +
theme_light()slice_max():slice_max() seleciona as n linhas maiores de acordo com uma das colunas.
slice_max() número de itens:Selecionando os países mais ricos em 2007 com o parâmetro n.
Depois aprenderemos como ordenar essas barras no gráfico
gapminder %>%
filter(year == 2007) %>%
slice_max(order_by = gdpPercap, n = 5) %>%
ggplot() +
geom_col(aes(x = country, y = gdpPercap)) +
theme_light() slice_max(), proporção das linhas:O parâmetro prop permite a seleção de uma proporção das linhas
ricos <- gapminder %>%
filter(year == 2007) %>%
slice_max(prop = .2, order_by = gdpPercap ) %>%
mutate(categoria = "Rico")
pobres <- gapminder %>%
filter(year == 2007) %>%
slice_min(prop = .2, order_by = gdpPercap) %>%
mutate(categoria = "Pobre")
bind_rows(ricos, pobres) %>%
ggplot(aes(x = gdpPercap, y = lifeExp)) +
geom_point( aes(color = continent)) +
facet_grid(. ~ categoria, scales = "free_x") +
geom_smooth(method = "lm", se = FALSE) +
theme_light() +
theme(
legend.position = "top"
)## `geom_smooth()` using formula 'y ~ x'
slice():slice() seleciona “linhas” de um data frame/tibble
slice():classificacao_brasileirao <- read_html("https://pt.wikipedia.org/wiki/Campeonato_Brasileiro_de_Futebol_de_2019_-_S%C3%A9rie_A") %>%
html_nodes("table") %>%
extract2(7) %>%
html_table()
limbo <- classificacao_brasileirao %>%
slice(12:16) %>%
select(time = Equipes )
limbo## time
## 1 Vasco da Gama
## 2 Atlético Mineiro
## 3 Fluminense
## 4 Botafogo
## 5 Ceará
group_by():group_by():A função group_by() será bastante usada.
Quem conhece SQL pode estranhar um pouco o comportamento desta função, pois ela não agrupa os dados diminuindo o número de linhas imediatamente.
Mas veja que ela indica que há agrupamento
Ela particiona o tibble. As operações passam a ser executadas em cada partição.
## # A tibble: 1,649 x 3
## # Groups: country [165]
## country date SI.POV.GINI
## <chr> <dbl> <dbl>
## 1 Albania 2017 33.2
## 2 Albania 2016 33.7
## 3 Albania 2015 32.9
## 4 Albania 2014 34.6
## 5 Albania 2012 29
## 6 Albania 2008 30
## 7 Albania 2005 30.6
## 8 Albania 2002 31.7
## 9 Albania 1996 27
## 10 Algeria 2011 27.6
## # ... with 1,639 more rows
group_by() (cont.):Para várias operações, o agrupamento faz com que o comportamento seja diferente.
Uma operação bastante usada é numerar as linhas de um tibble.
No tibble agrupado, essa operação acontece em cada grupo.
group_by() (cont.):As funções lag() no contexto da group_by e lead() no contexto da group_by funcionam dentro de cada grupo (o primeiro value de um grupo não acessa o valor do outro grupo com lag() .
summarise() (cont.):group_by() (cont.):A função group_by só leva a uma sumarização, ou seja, só transforma o tibble em um tibble com o número de linhas igual ao número de grupos, quando executamos a função summarise()
Se executarmos summarise() sem particionar o tibble, a operação resulta em uma linha.
maiores_temp_dia <- dados_heathrow %>%
group_by(date(date)) %>%
summarise(
maxima = max(air_temp),
minima = min(air_temp),
media = mean(air_temp)
)## `summarise()` ungrouping output (override with `.groups` argument)
A função slice_max() no contexto da group_by() retorna os n maiores valores. Se o tibble estiver agrupado, pra cada grupo.
maiores_temp_dia <- dados_heathrow %>%
group_by(date(date)) %>%
slice_max(n = 1, order_by = air_temp) %>%
ungroup() %>%
mutate(
hora = hour(date),
estacao =
case_when(
month(date) %in% 1:3 ~ "Inverno",
month(date) %in% 7:9 ~ "Verão",
TRUE ~ "Outono/Primavera"
)
) %>%
select(hora, estacao, air_temp)
ggplot(maiores_temp_dia) +
geom_density( aes(x = hora, color = estacao )) +
theme_light()A biblioteca dplyr chegou à versão 1.0.0 em junho de 2020, estabelecendo sua interface permanente e ganhando novas funcionalidades, como as column-wise operations, com a função across().
A função across(), quando usada em conjunto com mutate(), summarise(), filter()… , facilita a execução de operações iguais em diversos campos com base no nome dos campos ou no valor deles.
A função across() tem dois parâmetros.
.cols define em quais colunas serão aplicados as funções.fns define a função (ou a lista de funções) que serão aplicadas em cada colunacandidatos %>%
mutate(
across(
.cols = starts_with("DATA") ,
.fns = dmy
)
) %>%
relocate(starts_with("DATA") ) %>%
str() ## tibble [29,146 x 58] (S3: tbl_df/tbl/data.frame)
## $ DATA_GERACAO : Date[1:29146], format: "2020-10-01" "2020-10-01" ...
## $ DATA_ELEICAO : Date[1:29146], format: "2018-10-07" "2018-10-07" ...
## $ DATA_NASCIMENTO : Date[1:29146], format: "1974-04-08" "1994-11-11" ...
## $ HORA_GERACAO : 'hms' num [1:29146] 11:48:48 11:48:48 11:48:48 11:48:48 ...
## ..- attr(*, "units")= chr "secs"
## $ ANO_ELEICAO : num [1:29146] 2018 2018 2018 2018 2018 ...
## $ COD_TIPO_ELEICAO : num [1:29146] 2 2 2 2 2 2 2 2 2 2 ...
## $ NOME_TIPO_ELEICAO : chr [1:29146] "ELEIÇÃO ORDINÁRIA" "ELEIÇÃO ORDINÁRIA" "ELEIÇÃO ORDINÁRIA" "ELEIÇÃO ORDINÁRIA" ...
## $ NUM_TURNO : num [1:29146] 1 1 1 1 1 1 1 1 1 1 ...
## $ COD_ELEICAO : num [1:29146] 297 297 297 297 297 297 297 297 297 297 ...
## $ DESCRICAO_ELEICAO : chr [1:29146] "Eleições Gerais Estaduais 2018" "Eleições Gerais Estaduais 2018" "Eleições Gerais Estaduais 2018" "Eleições Gerais Estaduais 2018" ...
## $ ABRANGENCIA : chr [1:29146] "ESTADUAL" "ESTADUAL" "ESTADUAL" "ESTADUAL" ...
## $ SIGLA_UF : chr [1:29146] "AC" "AC" "AC" "AC" ...
## $ SIGLA_UE : chr [1:29146] "AC" "AC" "AC" "AC" ...
## $ DESCRICAO_UE : chr [1:29146] "ACRE" "ACRE" "ACRE" "ACRE" ...
## $ CODIGO_CARGO : num [1:29146] 7 6 7 6 6 7 7 7 7 7 ...
## $ DESCRICAO_CARGO : chr [1:29146] "DEPUTADO ESTADUAL" "DEPUTADO FEDERAL" "DEPUTADO ESTADUAL" "DEPUTADO FEDERAL" ...
## $ SEQUENCIAL_CANDIDATO : num [1:29146] 1e+10 1e+10 1e+10 1e+10 1e+10 ...
## $ NUMERO_CANDIDATO : num [1:29146] 11111 3131 12300 5555 2722 ...
## $ NOME_CANDIDATO : chr [1:29146] "GEHLEN DINIZ ANDRADE" "JEFFERSON BARROSO DE ARAÚJO " "ANTONIO FERREIRA PEREIRA" "MARIVALDO GONÇALVES DE MELO" ...
## $ NOME_URNA_CANDIDATO : chr [1:29146] "GERLEN DINIZ" "JEFFERSON BARROSO " "TONY FERREIRA - CANTOR" "MARIVALDO MELO" ...
## $ NOME_SOCIAL_CANDIDATO : chr [1:29146] "#NULO#" "#NULO#" "#NULO#" "#NULO#" ...
## $ CPF_CANDIDATO : chr [1:29146] "35954590249" "85546976268" "30843766204" "27608417234" ...
## $ EMAIL_CANDIDATO : chr [1:29146] "ANGELMFERREIRA@GMAIL.COM" "JEFFERSON.BARROSO.PROJETOS@GMAIL.COM" "TONY-FPEREIRA@HOTMAIL.COM" "MELOSENA@HOTMAIL.COM" ...
## $ COD_SITUACAO_CANDIDATURA : num [1:29146] 12 12 12 12 12 12 12 12 12 12 ...
## $ DES_SITUACAO_CANDIDATURA : chr [1:29146] "APTO" "APTO" "APTO" "APTO" ...
## $ COD_DETALHE_SITUACAO_CAND : num [1:29146] 2 2 2 2 2 2 2 2 2 2 ...
## $ DES_DETALHE_SITUACAO_CAND : chr [1:29146] "DEFERIDO" "DEFERIDO" "DEFERIDO" "DEFERIDO" ...
## $ TIPO_AGREMIACAO : chr [1:29146] "COLIGAÇÃO" "COLIGAÇÃO" "PARTIDO ISOLADO" "COLIGAÇÃO" ...
## $ NUMERO_PARTIDO : num [1:29146] 11 31 12 55 27 17 17 19 23 77 ...
## $ SIGLA_PARTIDO : chr [1:29146] "PP" "PHS" "PDT" "PSD" ...
## $ NOME_PARTIDO : chr [1:29146] "PROGRESSISTAS" "PARTIDO HUMANISTA DA SOLIDARIEDADE" "PARTIDO DEMOCRÁTICO TRABALHISTA" "PARTIDO SOCIAL DEMOCRÁTICO" ...
## $ CODIGO_LEGENDA : num [1:29146] 1e+10 1e+10 1e+10 1e+10 1e+10 ...
## $ NOME_COLIGACAO : chr [1:29146] "MUDANÇA COM COMPROMISSO" "Frente Popular do Acre I" "PARTIDO ISOLADO" "Mudança e Competência 1" ...
## $ COMPOSICAO_LEGENDA : chr [1:29146] "PP / PMN / PPS / PTC / PR" "PT / PSB / PHS / PC do B / DC" "PDT" "PP / PSDB / PSD / MDB / DEM / SOLIDARIEDADE / PTC / PMN / PR / PTB / PPS" ...
## $ CODIGO_NACIONALIDADE : num [1:29146] 1 1 1 1 1 1 1 1 1 1 ...
## $ DESCRICAO_NACIONALIDADE : chr [1:29146] "BRASILEIRA NATA" "BRASILEIRA NATA" "BRASILEIRA NATA" "BRASILEIRA NATA" ...
## $ SIGLA_UF_NASCIMENTO : chr [1:29146] "AC" "AC" "AC" "AC" ...
## $ CODIGO_MUNICIPIO_NASCIMENTO : num [1:29146] -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 ...
## $ NOME_MUNICIPIO_NASCIMENTO : chr [1:29146] "RIO BRANCO" "RIO BRANCO " "FEIJÓ" "RIO BRANCO" ...
## $ IDADE_DATA_POSSE : num [1:29146] 46 26 52 56 41 61 37 48 44 49 ...
## $ NUM_TITULO_ELEITORAL_CANDIDATO: chr [1:29146] "002092722470" "006267492410" "000385762402" "004960932291" ...
## $ CODIGO_SEXO : num [1:29146] 2 2 2 2 4 2 2 4 2 4 ...
## $ DESCRICAO_SEXO : chr [1:29146] "MASCULINO" "MASCULINO" "MASCULINO" "MASCULINO" ...
## $ COD_GRAU_INSTRUCAO : num [1:29146] 8 8 6 8 6 2 8 4 6 6 ...
## $ DESCRICAO_GRAU_INSTRUCAO : chr [1:29146] "SUPERIOR COMPLETO" "SUPERIOR COMPLETO" "ENSINO MÉDIO COMPLETO" "SUPERIOR COMPLETO" ...
## $ CODIGO_ESTADO_CIVIL : num [1:29146] 3 1 1 9 3 3 3 3 1 3 ...
## $ DESCRICAO_ESTADO_CIVIL : chr [1:29146] "CASADO(A)" "SOLTEIRO(A)" "SOLTEIRO(A)" "DIVORCIADO(A)" ...
## $ CODIGO_COR_RACA : chr [1:29146] "01" "01" "03" "01" ...
## $ DESCRICAO_COR_RACA : chr [1:29146] "BRANCA" "BRANCA" "PARDA" "BRANCA" ...
## $ CODIGO_OCUPACAO : num [1:29146] 277 999 257 257 257 999 999 999 999 160 ...
## $ DESCRICAO_OCUPACAO : chr [1:29146] "DEPUTADO" "OUTROS" "EMPRESÁRIO" "EMPRESÁRIO" ...
## $ DESPESA_MAX_CAMPANHA : num [1:29146] 1000000 2500000 1000000 2500000 2500000 1000000 1000000 1000000 1000000 1000000 ...
## $ COD_SIT_TOT_TURNO : num [1:29146] 2 5 5 5 5 5 5 5 5 5 ...
## $ DESC_SIT_TOT_TURNO : chr [1:29146] "ELEITO POR QP" "SUPLENTE" "SUPLENTE" "SUPLENTE" ...
## $ SITUACAO_REELEICAO : chr [1:29146] "S" "N" "N" "N" ...
## $ SITUACAO_DECLARAR_BENS : chr [1:29146] "S" "N" "S" "S" ...
## $ NUMERO_PROTOCOLO_CANDIDATURA : chr [1:29146] "-1" "-1" "-1" "-1" ...
## $ NUMERO_PROCESSO : chr [1:29146] "06006019120186010000" "06005594220186010000" "06003446620186010000" "06004320720186010000" ...
## - attr(*, ".internal.selfref")=<externalptr>
É possível usar a função across() em conjunto com filter. O código abaixo seleciona apenas as linhas pras quais nenhum percentil é nulo.
percentis_pisas <- wb_search(pattern = "LO\\.PISA\\.SCI\\.P.*")
dados_pisa <- wb_data(indicator = percentis_pisas$indicator_id )
pisas_completos <- dados_pisa %>%
filter(
across(starts_with("LO"), ~!is.na(.x) )
) %>%
left_join(
country_codes %>% rename(country_gapminder = country),
by = c("iso3c" = "iso_alpha")
) %>%
left_join(
gapminder %>% select(country, continent) %>% distinct(),
by = c("country")
) %>%
relocate(
continent
) %>%
rename_with(
~str_remove(.x, "LO\\.PISA\\.SCI\\.")
)
head(pisas_completos)## # A tibble: 6 x 14
## continent iso2c iso3c country date P05 P10 P25 P50 P75 P90 P95
## <fct> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Europe AL ALB Albania 2000 221 256 315 376. 438 497 531
## 2 Europe AL ALB Albania 2009 242. 276. 331. 391. 454. 504. 532.
## 3 Europe AL ALB Albania 2012 221. 271. 340. 397. 464. 517. 549.
## 4 Europe AL ALB Albania 2015 301. 328. 373. 426. 481. 530. 558.
## 5 Europe AL ALB Albania 2018 298. 323. 366. 416. 466. 514. 541.
## 6 <NA> AE ARE United~ 2012 299. 328. 382. 448. 512. 572. 605.
## # ... with 2 more variables: country_gapminder <chr>, iso_num <int>
Podemos executar um summarise across, selecionando as colunas desejadas
Neste caso, tiramos, para o último ano disponível, a mediana dos percentis para cada continente
sumario <- pisas_completos %>%
select(
c(continent, date, country) | matches("^P[0-9]{2}$")
) %>%
filter(
date == max(date),
!is.na(continent)
) %>%
select(
-date
) %>%
group_by(continent) %>%
summarise(
across(
.cols = matches("P[0-9]{2}"),
.fns = median
)
) %>%
view()## `summarise()` ungrouping output (override with `.groups` argument)
O parâmetro .fns aceita uma lista nomeada de funções.
É possível executar várias operações de uma só vez criando novos campos. O nome do elemento relativo à função na lista é concatenado com os campos originais gerando um novo nome pros resultados das operações
sumario <- pisas_completos %>%
select(
c(continent, date, country) | matches("P[0-9]{2}")
) %>%
filter(
!is.na(continent)
) %>%
summarise(
across(
.cols = matches("P[0-9]{2}"),
.fns = list(mediana = median, media = mean)
)
)
sumario %>%
gt() %>%
fmt_number(
columns = matches("P[0-9]{2}"),
decimals = 2
)| P05_mediana | P05_media | P10_mediana | P10_media | P25_mediana | P25_media | P50_mediana | P50_media | P75_mediana | P75_media | P90_mediana | P90_media | P95_mediana | P95_media |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 327.10 | 320.71 | 362.56 | 353.52 | 423.96 | 410.09 | 491.12 | 473.15 | 557.80 | 537.27 | 610.70 | 590.77 | 641.69 | 621.33 |
Citamos anteriormente que o tibble/dataframe pode conter colunas que são listas de elementos.
É muito útil armazenar dataframes aninhados nestas list-columns, ou seja, uma coluna que contém um dataframe por linha.
Esta técnica tem inúmeros usos. Por exemplo, rodar modelos em partições dos dados originais
No código seguinte, particionamos os dados do PISA por país e aninhamos em uma list-column usando a função nest(). Podemos dar group_by e depois nest ou explicitar as colunas que queremos aninhar.
É uma boa prática rodar ungroup assim que não for necessário usar o tibble agrupado
pisa_nested <- pisas_completos %>%
select(country, date, P05) %>%
arrange(date) %>%
group_by(country) %>%
mutate(
ultimo_P05 = last(P05)
) %>%
ungroup() %>%
nest(pisa = -c(country, ultimo_P05)) %>%
ungroup()
head(pisa_nested)## # A tibble: 6 x 3
## country ultimo_P05 pisa
## <chr> <dbl> <list>
## 1 Albania 298. <tibble [5 x 2]>
## 2 Argentina 261. <tibble [6 x 2]>
## 3 Australia 334. <tibble [7 x 2]>
## 4 Belgium 328. <tibble [7 x 2]>
## 5 Bulgaria 279. <tibble [6 x 2]>
## 6 Brazil 268. <tibble [7 x 2]>
Agora podemos rodar uma regressão linear, por exemplo.
Vamos usar programação funcional para rodar lm(), a função de regressão linear, para cada país
Temos, agora, o resultado da execução de um modelo de regressão linear pra cada país.
Precisamos, pra cada objeto retornado, extrair os coeficientes.
vamos usar a função broom::tidy()
Essa função não está preparada para receber um vetor de objetos mas sim UM único objeto.
Portanto, essa operação deve ser feita linha a linha. Para que isso aconteça é necessário rodar a função rowwise, que equivale a um group_by() por linha.
A função broom::tidy() criou novos tibbles aninhados com as informações a respeito dos coeficientes estimados.
Vamos, agora explodir esses tibbles, com a função unnest(), mostrando o coeficiente que nos interessa aqui, o \(\beta\) relativo à variável explicativa contendo o ano.
pisa_modelo_coefs_exp <- pisa_modelo_coefs %>%
unnest(coeficientes) %>%
filter(
term == "date"
) %>%
select(
country,
ultimo_P05,
estimate,
std.error,
p.value
) %>%
filter(
!is.na(estimate)
) %>%
arrange(
estimate %>% desc()
) %>%
mutate(
posicao = row_number(),
.before = everything()
) %>%
slice(
c(1:10,(nrow(.)-10+1):nrow(.))
)
gt(pisa_modelo_coefs_exp) %>%
fmt_number(
column = "ultimo_P05",
decimals = 0
) %>%
fmt_number(
column = c("estimate","std.error", "p.value"),
decimals = 2
)| posicao | country | ultimo_P05 | estimate | std.error | p.value |
|---|---|---|---|---|---|
| 1 | Trinidad and Tobago | 279 | 7.51 | NaN | NaN |
| 2 | Peru | 280 | 5.52 | 0.64 | 0.00 |
| 3 | Liechtenstein | 383 | 5.36 | 0.95 | 0.01 |
| 4 | Argentina | 261 | 4.62 | 1.53 | 0.04 |
| 5 | Albania | 298 | 4.45 | 2.15 | 0.13 |
| 6 | Vietnam | 418 | 3.41 | 0.76 | 0.14 |
| 7 | Malaysia | 313 | 3.41 | 3.20 | 0.48 |
| 8 | Colombia | 287 | 3.40 | 0.78 | 0.02 |
| 9 | Kazakhstan | 284 | 3.35 | 5.54 | 0.61 |
| 10 | Qatar | 259 | 3.34 | 1.66 | 0.14 |
| 66 | Netherlands | 329 | −2.25 | 0.55 | 0.02 |
| 67 | Finland | 356 | −2.33 | 1.04 | 0.08 |
| 68 | Croatia | 327 | −2.35 | 0.58 | 0.03 |
| 69 | Slovak Republic | 307 | −2.57 | 0.97 | 0.06 |
| 70 | Kyrgyz Republic | 183 | −2.58 | NaN | NaN |
| 71 | Costa Rica | 300 | −2.59 | 0.41 | 0.10 |
| 72 | Georgia | 255 | −3.98 | NaN | NaN |
| 73 | Lebanon | 237 | −4.04 | NaN | NaN |
| 74 | United Arab Emirates | 272 | −4.51 | 0.38 | 0.05 |
| 75 | Azerbaijan | 257 | −14.07 | NaN | NaN |
Até agora acessamos dados que estavam disponíveis em bibliotecas, mas muitas vezes encontramos dados em arquivos.
De modo geral, as funções da biblioteca readr são mais rápidas do que as da biblioteca base, e além disso mostram barra de progresso no console. É possível reconhecê-las pelo _ ao invés de .
O caso mais comum é ler dados em formato de tabela para um tibble
Fonte: (RStudio 2019a)
O portal da CVM é uma fonte interessante de dados
O código abaixo baixa os dados que ainda não estão na nossa base. E, por via das dúvidas, baixa o último arquivo de novo.
Note que vamos usar programação funcional no exemplo. Sem uso de loop vamos executar uma função para cada linha do data frame/tibble.
Esta função vai salvar os dados como um objeto rds. Este formato é binário, diferentemente do formato csv (textual). Podemos compactar estes dados ou não. Se usarmos compactação, a leitura e escrita são mais lentas, mas os arquivos ocupam menos espaço.
existem <- tibble(arquivo = list.files("dados/fundos")) %>%
slice_head(n = nrow(.)-1)
salva <- function(endereco, arquivo){
print(endereco)
read_csv2(endereco) %>%
mutate(
across(
-c("CNPJ_FUNDO", "DT_COMPTC"),
as.numeric
)
) %>%
write_rds(paste0("dados/fundos/",arquivo), compress = "xz")
}
baixa_faltantes <- tibble(data_dado = seq(ymd("2019-01-01"), by = "month", ymd(today()-4) )) %>%
mutate(data_formato = format(data_dado, "%Y%m")) %>%
mutate(
endereco = paste0(
"http://dados.cvm.gov.br/dados/FI/DOC/INF_DIARIO/DADOS/",
"inf_diario_fi_",
data_formato,
".csv")) %>%
mutate(arquivo = paste0("inf_diario_fi_",data_formato,".rds")) %>%
anti_join(existem, by = c("arquivo" = "arquivo")) %>%
select(-data_formato) %>%
mutate(data = pmap(.l = list(endereco = endereco, arquivo = arquivo), salva))## [1] "http://dados.cvm.gov.br/dados/FI/DOC/INF_DIARIO/DADOS/inf_diario_fi_202101.csv"
Após a operação de baixar os dados dos arquivos, podemos ler todos os arquivos de uma pasta e colocá-los em um só tibble, já que eles são do mesmo formato.
Usamos execução em paralelo no exemplo, que é possibilitada ao chamarmos uma função future_
Após a leitura dos arquivos eles ficam armazenados de forma aninhada numa coluna chamada “conteudo.” “conteudo” é uma coluna feita de tibbles, cada tibble correspondendo a uma linha do tibble principal. A coluna “conteudo” é uma nested list-column.
plan(multiprocess)
todos_os_fundos <- tibble(arquivo = list.files("dados/fundos")) %>%
mutate(arquivo = str_glue("dados/fundos/{arquivo}")) %>%
mutate(conteudo = future_map(arquivo, read_rds , .progress = TRUE)) %>%
unnest(conteudo)
head(todos_os_fundos)## # A tibble: 6 x 9
## arquivo CNPJ_FUNDO DT_COMPTC VL_TOTAL VL_QUOTA VL_PATRIM_LIQ CAPTC_DIA
## <glue> <chr> <date> <dbl> <dbl> <dbl> <dbl>
## 1 dados/~ 00.017.02~ 2019-01-02 1.13e8 2.67e13 113068421 0
## 2 dados/~ 00.017.02~ 2019-01-03 1.13e8 2.67e13 113078541 0
## 3 dados/~ 00.017.02~ 2019-01-04 1.14e8 2.67e13 113088700 0
## 4 dados/~ 00.017.02~ 2019-01-07 1.14e8 2.67e13 113098833 0
## 5 dados/~ 00.017.02~ 2019-01-08 1.14e8 2.67e13 113108977 0
## 6 dados/~ 00.017.02~ 2019-01-09 1.14e8 2.67e13 113119100 0
## # ... with 2 more variables: RESG_DIA <dbl>, NR_COTST <dbl>
A função unnest pode ser usada, como já vimos
Existem duas dois tipos de atividades envolvidas na leitura de páginas web:
Extração do conteúdo que interessa a partir da página já obtida
Automatização da navegação com uso de emuladores de browser
Vamos ver primeiro como ler o conteúdo de uma página e depois como usar um emulador para navegar em um portal e obter o que desejamos
Para efetuar a leitura de páginas WEB, é necessário conhecer como é a estrutura de uma página web.
Uma página web pode ser representada por uma árvore de objetos, também chamada de DOM (Document Object Model).
Esta árvore de objetos é definida pelo conteúdo de linguagem html que existe na página (e pode ser modificado por scripts em javascript e definições de estilo do CSS).
Podem existir objetos de vários tipos em uma página: links, inputs de dados, tabelas, células de tabelas, parágrafos, cabeçalhos etc.
Para retirar da página web o conteúdo de que precisamos, temos que analisar como é esta árvore de objetos e que nós desta árvore nos interessam.
Imagine que queremos buscar dados na página de histórico de preços e taxas dos títulos brasileiros
A tecla F12 do Chrome nos permite ver a árvore DOM da página em que estamos navegando.
É possível clicar com o botão direito e inspecionar um elemento específico de forma a saber onde ele está na árvore e que tipo de elemento ele é (mesmo que você saiba pouco de html).
O mais importante é saber que uma tag html que define um elemento tem a sintaxe:
<tipo_elemento nome_atributo=valor_attributo>texto do elemento</tipo_elemento>
Mesmo sem saber hmtl, fica claro que queremos esse tal de atributo href dos tais elementos a seja lá o que diabos isso seja (a é um link e href é o destino do link).
A biblioteca rvest possibilita a extração destes elementos.
É possível caminhar pela árvore DOM até os nós desejados e atributos que queremos usando html_nodes() e html_attr().
Munidos de uma função que faz download e salva um arquivo, podemos caminhar pelas planilhas e salvá-las
existem <- tibble(arquivo = list.files("dados/titulos"))
salva_planilha <- function(name, endereco){
arquivo <- endereco
destino <- name
print(arquivo)
download.file(
arquivo,
paste0("dados/titulos/",destino ),
mode = "wb" ##PRA ARQUIVOS BINÁRIOS,
)
}
links <- read_html("https://sisweb.tesouro.gov.br/apex/f?p=2031:2:0::::") %>%
html_nodes("body") %>%
html_nodes("a") %>%
html_attr("href") %>%
enframe(value = "endereco") %>%
filter(str_detect(endereco, "cosis/sistd/obtem_arquivo/")) %>%
mutate(destino = str_extract(endereco, "[0-9]*:[0-9]*$") %>% str_remove(":") ) %>%
mutate(destino = paste0(destino,".xls")) %>%
anti_join(existem, by = c("destino" = "arquivo")) %>%
mutate(endereco = paste0("https://sisweb.tesouro.gov.br/apex/",endereco) ) %>%
mutate(data = map2(destino, endereco, salva_planilha ))vamos aproveitar as planilhas que baixamos para ver como podemos efetuar a leitura de planilhas Excel.
Agora vamos organizar as planilhas que lemos do site do tesouro.
Vamos usar a biblioteca read_xl.
Precisamos ler todas as sheets de todos os arquivos em um diretório (onde baixamos os arquivos excel do site do tesouro).
A função abaixo lê as sheets de um arquivo.
A função abaixo lê o conteúdo de cada sheet
Munidos destas duas funções, podemos guardar tudo em um só data frame.
Primeiro, definindo as sheets a ler
plan(multiprocess)
sheets_pra_ler <- list.files("dados/titulos") %>%
enframe(value = "arquivo") %>%
mutate(sheet = future_map(arquivo, le_sheets, .progress = TRUE)) %>%
unnest(cols = sheet) Depois lendo as sheets para um data frame
Neste momento, vamos executar a função em processamento paralelo, usando todos os processadores da nossa máquina.
Para isso, vamos lançar mão da biblioteca furrr%.
Ela contém adaptações de todas as funções map. São as funções baseadas na future_map.
Antes de rodar a future_map, determinamos a forma de paralelização. No nosso caso, plan(multiprocess).
É legal notar a forma como usamos expressão regular aqui
plan(multiprocess)
sheets_lidas <- sheets_pra_ler %>%
mutate(titulo = str_extract(sheet,"[^[0-9]]*" )) %>%
mutate(vencimento = str_extract(sheet,"[0-9]{6}" )) %>%
mutate(vencimento = dmy(vencimento)) %>%
mutate(
arquivo_out = arquivo,
sheet_out = sheet
) %>%
mutate(data = future_map2(.x = arquivo, .y = sheet , le_conteudo_sheet, .progress = TRUE)) %>%
unnest(data) ##
Progress: --------------------------------------------------------------------------------------------------------------- 100%
Uma última arrumada
taxas_titulos <- sheets_lidas %>%
rename(
data = 8,
taxa_compra_manha = 9,
taxa_venda_manha = 10,
pu_compra_manha = 11,
pu_venda_manha = 12,
pu_base_manha = 13
) %>%
mutate(
data = if_else(
str_detect(data, "/"),
dmy(data),
as.numeric(data) + ymd("1899-12-31")
)
) %>%
mutate(
titulo = str_trim(titulo),
taxa_compra_manha = as.numeric(taxa_compra_manha),
taxa_venda_manha = as.numeric(taxa_venda_manha),
pu_compra_manha = as.numeric(pu_compra_manha),
pu_venda_manha = as.numeric(pu_venda_manha),
pu_base_manha = as.numeric(pu_base_manha)
) %>%
select(
titulo,
vencimento,
data,
taxa_compra_manha,
taxa_venda_manha,
pu_compra_manha,
pu_venda_manha,
pu_base_manha
) %>%
mutate(
titulo = if_else(
titulo == "NTN-B Principal",
"NTN-B Princ",
titulo
)
)Aí fica fácil fazer a análise que desejarmos
ntnb_2045 <- taxas_titulos %>%
filter(
titulo == "NTN-B Princ",
vencimento == ymd("2045-05-15")
) %>%
mutate(taxa = (taxa_compra_manha + taxa_venda_manha)/2 )
ggplot(ntnb_2045) +
geom_line(aes(x = data, y = taxa) ) +
scale_x_date(
date_breaks = "3 months",
limits = c(ymd("2016-12-01"),NA),
labels = date_format("%m/%y")
) +
scale_y_continuous(
labels = percent_format(accuracy = 0.1)
) +
labs(y = "NTN-B Principal 2045", x = "Data") +
theme_light()Neste caso, precisamos entender quando requisição é feita pela página e emular essas requisições passando os parâmetros que devem ser passados via método POST.
Primeiro usamos a função crossing para criar um novo tibble que combina todas as \(n\) linhas do tibble com os grupos e todas as \(m\) linhas do tibble de datas, resultando num tibble com \(m * n\) linhas, pois queremos pegar o resultado de todos grupos para todas as datas.
url <- "http://www.ceagesp.gov.br/entrepostos/servicos/cotacoes/#cotacao"
le_pagina_ceagesp <- function(grupo, data){
print(str_glue("Grupo:{grupo}, data:{data}"))
dados <- NA
fd <- list(
cot_grupo = grupo,
cot_data = data
)
resp <- POST(url, body=fd, encode="form")
tabela <- content(resp) %>%
html_nodes("table") %>%
extract2(1) %>%
html_table()
nomes <- tabela %>% slice(2)
dados <- tabela %>% slice(3:nrow(tabela)) %>%
mutate(data = dmy(data))
names(dados) <- c(nomes, "data_precos")
dados
}
le_pagina_ceagesp_ignora_erro <- possibly(le_pagina_ceagesp, otherwise = NA)
grupos <- c(
"FRUTAS",
"LEGUMES",
"VERDURAS",
"DIVERSOS",
"FLORES",
"PESCADOS"
)
datas <- seq(from = today()-1, by = 1, to = today()-1) %>%
format("%d/%m/%Y")
dados_ceagesp <- enframe(grupos, value = "grupo", name = "id_grupo") %>%
crossing(enframe(datas, value = "data", name = "id_data")) %>%
mutate(leitura = map2(.x = grupo, .y = data, le_pagina_ceagesp_ignora_erro )) %>%
filter(!is.na(leitura)) %>%
unnest(cols = leitura)## Grupo:FRUTAS, data:01/12/2020
## Grupo:LEGUMES, data:01/12/2020
## Grupo:VERDURAS, data:01/12/2020
## Grupo:DIVERSOS, data:01/12/2020
## Grupo:FLORES, data:01/12/2020
## Grupo:PESCADOS, data:01/12/2020
## # A tibble: 0 x 5
## # ... with 5 variables: id_grupo <int>, grupo <chr>, id_data <int>, data <chr>,
## # leitura <???>
Principalmente em dois casos precisamos emular um browser para obter o conteúdo html que desejamos, para então usar a rvest se necessário:
Quando a página precisa executar um script javascript para que fique pronta
Quando a página exige muitos passos de navegação antes de chegarmos aos dados desejados;
É preciso estudar cada caso para entender qual melhor estratégia. Precisamos pensar em
Estabilidade: quão dependente a leitura é de aspectos específicos dos dados atuais ou da forma como ela foi desenvolvida
Correção: temos que nos certificar de que o resultado que obtemos é o mesmo que a página está disponibilizando ao usuário. A página muitas vezes é criada para atender um usuário humano. Uma página que é carregada sem dados ou com dados preliminares que depois são preeenchidos ou alterados por um script é adequada ara humanos mas pode não ser adequada para uma extração automática.
Selenium é um famoso automatizador de browser. Através dele, é possível enviar cliques, preencher caixas de texto etc. até se chegar aos dados desejados. Além disso, ele executa os scripts javascripts que a página exige.
A biblioteca RSelenium possibilita o uso na linguagem R
Abaixo um exemplo de ativação do firefox a partir do R com RSelenium
library(RSelenium)
library(rvest)
library(httr)
library(magrittr)
# initiate selenium -- start it this way every time
# make sure the location of firefox on your machine is correct
rs <- rsDriver(
browser = "firefox",
extraCapabilities = list(
`mox:firefoxOptions` = list(
binary = "C:/Program Files/Mozilla Firefox/firefox.exe"
)
)
)
preparou <- function(){
outputzipfile <- "nada"
try(outputzipfile <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[2]/div[2]/div/div[1]"))
if(typeof(outputzipfile) == "S4"){
resposta <- outputzipfile$getElementText() == "Output Zip File"
}
else{
resposta <- FALSE
}
resposta
}
rsc <- rs$client
rsc$navigate("https://sigel.aneel.gov.br/Down/")
Sys.sleep(3)
imagem <- rsc$findElement(using = "css selector", "#uniqName_9_1 > div:nth-child(1)" )
imagem$clickElement()
Sys.sleep(3)
UHE <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[1]/div[1]/div[1]/div[2]/div/div[3]/div[1]")
UHE$clickElement()
GD <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[1]/div[1]/div[1]/div[2]/div/div[1]/div[1]")
GD$clickElement()
baixar <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[1]/div[2]/div[2]/div")
baixar$clickElement()
while(!preparou()){
print("esperando")
Sys.sleep(1)
}
arquivo_download <- rsc$findElement(using = "xpath", "/html/body/div[2]/div[1]/div/div[13]/div[2]/div/div/div[3]/div[2]/div[2]/div[2]/div/div[2]/div/a")
endereco <- arquivo_download$getElementText()
download.file(url = unlist(endereco), destfile = "dados/aneel/gd.zip")O data frame mostrad neste slide não está “Tidy,” ou seja não está de forma que cada linha represente um evento e cada coluna represente um atributo do evento.
Para nós, neste caso, um evento é formado por todas as informações de um ativo em um dia e não uma só das informações de um ativo em um dia.
Isso porque é extremamente comum fazermos contas com mais de uma informação do ativo em um dia (máxima - mínima, por exemplo)
dados_long <- tribble(
~data, ~ativo, ~cotacao, ~valor,
make_date(2010,01,01), "PETR4", "Máxima", 1.02,
make_date(2010,01,01), "PETR4", "Mínima", 1.00,
make_date(2010,01,02), "PETR4", "Máxima", 1.05,
make_date(2010,01,02), "PETR4", "Mínima", 1.02,
make_date(2010,01,03), "PETR4", "Máxima", 1.07,
make_date(2010,01,03), "PETR4", "Mínima", 1.01,
make_date(2010,01,04), "PETR4", "Máxima", 1.07,
make_date(2010,01,04), "PETR4", "Mínima", 1.03
)O pacote Tidyr ajuda a arrumar os data frames dessas formas. O hex sticker dele é bem explicativo.
pivot_longer() e pivot_wider()Os nomes das principais funções mudaram em setembro de 2019 (quando saiu a versão 1.0.0). Antes se chamavam gather() e spread() e agora se chamam pivot_wider() e pivot_longer(), o que é mais intuitivo.
balanços
download.file(
url = "http://dados.cvm.gov.br/dados/CIA_ABERTA/DOC/DFP/DADOS/dfp_cia_aberta_2019.zip",
destfile = "dados/info_empresas_cvm/info.zip"
)
unzip(zipfile = "dados/info_empresas_cvm/info.zip", exdir = "dados/info_empresas_cvm/unzip" )A função que faz isso se chama pivot_wider()
Seus parâmetros mais usados são:
data, que é o data frame a ser tratado
names_from, que é o atributo de onde vêm os nomes para os novos atributos
values_from, o atributo de onde vêm os valores dos novos atributos
dados_wide <- dados_long %>%
pivot_wider(
names_from = cotacao,
values_from = valor
)
str(dados_wide)## tibble [4 x 4] (S3: tbl_df/tbl/data.frame)
## $ data : Date[1:4], format: "2010-01-01" "2010-01-02" ...
## $ ativo : chr [1:4] "PETR4" "PETR4" "PETR4" "PETR4"
## $ Máxima: num [1:4] 1.02 1.05 1.07 1.07
## $ Mínima: num [1:4] 1 1.02 1.01 1.03
É muito comum recebermos os dados com colunas que deviam ser valores de um atributo, e não um atributo em si.
O exemplo clássico é a colocação de datas nas colunas do dado, como nos dados retirados do site Datasus
read_csv2(
"dados/siab/cadastro_numero_familias.csv",
skip = 3,
locale = locale(encoding = "latin1" )
) %>%
glimpse()## Rows: 5,502
## Columns: 19
## $ Município <chr> "110001 Alta Floresta D'Oeste", "110037 Alto Alegre dos P...
## $ `1998` <chr> "2", "1224", "-", "3797", "676", "-", "-", "-", "5291", "...
## $ `1999` <chr> "3458", "2204", "546", "2763", "9922", "2268", "-", "628"...
## $ `2000` <chr> "4668", "1911", "2002", "2832", "11790", "2615", "1749", ...
## $ `2001` <chr> "5078", "1994", "2348", "3924", "11648", "3061", "1839", ...
## $ `2002` <chr> "5078", "1976", "2348", "3907", "11654", "3732", "1839", ...
## $ `2003` <chr> "5077", "1804", "2021", "3463", "9931", "5499", "1382", "...
## $ `2004` <chr> "4925", "1766", "932", "4077", "18830", "7580", "1803", "...
## $ `2005` <chr> "5865", "1698", "1956", "4141", "14531", "7470", "1859", ...
## $ `2006` <chr> "6518", "3623", "1994", "4244", "16509", "6459", "1941", ...
## $ `2007` <chr> "6489", "3500", "2922", "1642", "17094", "7502", "1929", ...
## $ `2008` <chr> "7093", "3682", "3195", "4493", "16936", "7712", "2381", ...
## $ `2009` <chr> "6946", "3172", "3210", "4693", "16948", "7546", "2344", ...
## $ `2010` <chr> "6836", "3292", "1480", "1586", "16528", "7112", "2592", ...
## $ `2011` <chr> "6918", "3416", "2551", "3122", "16912", "7378", "1855", ...
## $ `2012` <chr> "6699", "3448", "3175", "4162", "18109", "7385", "1962", ...
## $ `2013` <chr> "6336", "3617", "3029", "4242", "19270", "9730", "1945", ...
## $ `2014` <chr> "6201", "3343", "3438", "4242", "14661", "9084", "1939", ...
## $ `2015` <chr> "6181", "3280", "-", "-", "6307", "-", "-", "-", "19716",...
Acontece que “2009” não é um atributo, mas sim o valor de um atributo que deveria ser data
A função pivot_longer() faz a operação de que precisamos.
data, que é o data frame a ser tratado
names_from, que é o atributo de onde vêm os nomes para os novos atributos
values_from, o atributo de onde vêm os valores dos novos atributos
Note também a função separate(), que divide colunas de acordo com caracteres separadores.
siab_familias <- read_csv2(
"dados/siab/cadastro_numero_familias.csv",
skip = 3,
locale = locale(encoding = "latin1" )
) %>%
pivot_longer(
cols = -`Município`,
names_to = "data",
values_to = "familias"
) %>%
rename(
municipio = `Município`
) %>%
separate(
col = municipio,
into = c("cod_municipio", "municipio"),
sep = " ",
extra = "merge"
) %>%
mutate(
cod_municipio = as.integer(cod_municipio),
data = as.integer(data),
familias = as.integer(familias)
)
head(siab_familias)## # A tibble: 6 x 4
## cod_municipio municipio data familias
## <int> <chr> <int> <int>
## 1 110001 Alta Floresta D'Oeste 1998 2
## 2 110001 Alta Floresta D'Oeste 1999 3458
## 3 110001 Alta Floresta D'Oeste 2000 4668
## 4 110001 Alta Floresta D'Oeste 2001 5078
## 5 110001 Alta Floresta D'Oeste 2002 5078
## 6 110001 Alta Floresta D'Oeste 2003 5077
Para pegar mais informações dos municipios, vamos ler um arquivo baixado do IBGE
municipios <- read_excel("dados/ibge/populacao.xlsx", skip = 2, col_names = TRUE)
head(municipios, n = 30)## # A tibble: 30 x 4
## Cód. Município Ano `Tabela 6579 - População residente est~
## <chr> <chr> <chr> <dbl>
## 1 1100015 Alta Floresta D'Oeste ~ 2001 26919
## 2 <NA> <NA> 2002 27237
## 3 <NA> <NA> 2003 27563
## 4 <NA> <NA> 2004 29001
## 5 <NA> <NA> 2005 28629
## 6 <NA> <NA> 2006 29005
## 7 <NA> <NA> 2008 24577
## 8 <NA> <NA> 2009 24354
## 9 <NA> <NA> 2011 24228
## 10 <NA> <NA> 2012 24069
## # ... with 20 more rows
Ops… células mescladas
Não deixe o ódio tomar você…
A função fill() pode te ajudar
municipios_ajeitado <- municipios %>%
rename(
cod_municipio = 1,
nome_municipio = 2,
ano = 3,
populacao = 4
) %>%
fill(cod_municipio, .direction = "down") %>%
fill(nome_municipio, .direction = "down") %>%
mutate(UF = str_extract(nome_municipio,"\\([A-Z]{2}\\)")) %>%
mutate(UF = str_extract(UF,"[A-Z]{2}")) %>%
mutate(
cod_municipio = as.integer(cod_municipio),
ano = as.integer(ano)
)
head(municipios_ajeitado)## # A tibble: 6 x 5
## cod_municipio nome_municipio ano populacao UF
## <int> <chr> <int> <dbl> <chr>
## 1 1100015 Alta Floresta D'Oeste (RO) 2001 26919 RO
## 2 1100015 Alta Floresta D'Oeste (RO) 2002 27237 RO
## 3 1100015 Alta Floresta D'Oeste (RO) 2003 27563 RO
## 4 1100015 Alta Floresta D'Oeste (RO) 2004 29001 RO
## 5 1100015 Alta Floresta D'Oeste (RO) 2005 28629 RO
## 6 1100015 Alta Floresta D'Oeste (RO) 2006 29005 RO
Para juntar as informações de dois tibbles em um só, podemos fazer isso de três formas
Eis as funções de join de data frames
Fonte: (RStudio 2019b)
Anteriormente, pegamos informações do cadastro de famílias…
## # A tibble: 6 x 4
## cod_municipio municipio data familias
## <int> <chr> <int> <int>
## 1 110001 Alta Floresta D'Oeste 1998 2
## 2 110001 Alta Floresta D'Oeste 1999 3458
## 3 110001 Alta Floresta D'Oeste 2000 4668
## 4 110001 Alta Floresta D'Oeste 2001 5078
## 5 110001 Alta Floresta D'Oeste 2002 5078
## 6 110001 Alta Floresta D'Oeste 2003 5077
E de municípios
## # A tibble: 6 x 5
## cod_municipio nome_municipio ano populacao UF
## <int> <chr> <int> <dbl> <chr>
## 1 1100015 Alta Floresta D'Oeste (RO) 2001 26919 RO
## 2 1100015 Alta Floresta D'Oeste (RO) 2002 27237 RO
## 3 1100015 Alta Floresta D'Oeste (RO) 2003 27563 RO
## 4 1100015 Alta Floresta D'Oeste (RO) 2004 29001 RO
## 5 1100015 Alta Floresta D'Oeste (RO) 2005 28629 RO
## 6 1100015 Alta Floresta D'Oeste (RO) 2006 29005 RO
Os códigos são diferentes, mas
de_para_codigos <- read_csv2(
"http://blog.mds.gov.br/redesuas/wp-content/uploads/2018/06/Lista_Munic%C3%ADpios_com_IBGE_Brasil_Versao_CSV.csv",
locale = locale(encoding = "latin1")
) %>%
select(IBGE, IBGE7)
head(de_para_codigos)## # A tibble: 6 x 2
## IBGE IBGE7
## <dbl> <dbl>
## 1 110001 1100015
## 2 110002 1100023
## 3 110003 1100031
## 4 110004 1100049
## 5 110005 1100056
## 6 110006 1100064
Agora juntamos as informações do Siab com as informações dos municípios
siab_familias %>%
inner_join(de_para_codigos, by = c("cod_municipio" = "IBGE")) %>%
inner_join(municipios_ajeitado, by = c("IBGE7" = "cod_municipio", "data" = "ano") ) ## # A tibble: 71,500 x 8
## cod_municipio municipio data familias IBGE7 nome_municipio populacao UF
## <dbl> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 110001 Alta Flor~ 2001 5078 1.10e6 Alta Floresta~ 26919 RO
## 2 110001 Alta Flor~ 2002 5078 1.10e6 Alta Floresta~ 27237 RO
## 3 110001 Alta Flor~ 2003 5077 1.10e6 Alta Floresta~ 27563 RO
## 4 110001 Alta Flor~ 2004 4925 1.10e6 Alta Floresta~ 29001 RO
## 5 110001 Alta Flor~ 2005 5865 1.10e6 Alta Floresta~ 28629 RO
## 6 110001 Alta Flor~ 2006 6518 1.10e6 Alta Floresta~ 29005 RO
## 7 110001 Alta Flor~ 2008 7093 1.10e6 Alta Floresta~ 24577 RO
## 8 110001 Alta Flor~ 2009 6946 1.10e6 Alta Floresta~ 24354 RO
## 9 110001 Alta Flor~ 2011 6918 1.10e6 Alta Floresta~ 24228 RO
## 10 110001 Alta Flor~ 2012 6699 1.10e6 Alta Floresta~ 24069 RO
## # ... with 71,490 more rows
## # A tibble: 6 x 4
## cod_municipio municipio data familias
## <int> <chr> <int> <int>
## 1 110001 Alta Floresta D'Oeste 1998 2
## 2 110001 Alta Floresta D'Oeste 1999 3458
## 3 110001 Alta Floresta D'Oeste 2000 4668
## 4 110001 Alta Floresta D'Oeste 2001 5078
## 5 110001 Alta Floresta D'Oeste 2002 5078
## 6 110001 Alta Floresta D'Oeste 2003 5077
Seria legal saber qual o tamanho médio das famílias.
Este procedimento pode ficar como exercício, com o uso do dado abaixo
numero_medio_familias <- read_excel("dados/ibge/pessoas_por_familia.xlsx", skip = 2) %>%
select(c(1, 3, 7)) %>%
rename(
cod_municipio = 1,
pessoas = 2,
numero = 3
)
numero_medio_familias## # A tibble: 83,476 x 3
## cod_municipio pessoas numero
## <chr> <chr> <chr>
## 1 1100015 1 pessoa ...
## 2 <NA> 2 pessoas 1892
## 3 <NA> 3 pessoas 2008
## 4 <NA> 4 pessoas 1736
## 5 <NA> 5 pessoas 685
## 6 <NA> 6 pessoas 234
## 7 <NA> 7 pessoas 146
## 8 <NA> 8 pessoas 19
## 9 <NA> 9 pessoas 24
## 10 <NA> 10 pessoas 9
## # ... with 83,466 more rows
Nos próximos slides vemos antipadrões de visualização e dados.
Fonte: (Exame 2018)
Às vezes mesmo sem querer (será?)
PS. gráfico de 2014
Fonte: (Alves 2014)
Gráfico com sombra, 3D, estilos que dificultam a interpretação
Fonte: (Healy 2018)
O pessoal que usa Excel muitas vezes ama 3D, mas…
(Healy 2018)
Fonte: (White 2015)
Assim como os restaurantes devemos servir pizzas de dois sabores no máximo
Tentem descobrir qual a menor e a menor categoria em cada caso
Fonte: (Viz 2018)
Fonte: (Viz 2018)
Gráficos de dois eixos podem mostrar relações espúrias muito facilmente. E elas dependem da escolha da escala e dos limites dos eixos.
Fonte (Rost 2018)
O mesmo dado pode levar aos seguintes gráficos (e suas soluções)
Algumas dicas para boa visualização de dados:
Enfatize o dado relevante
Integre texto e gráfico
Use fontes e cores diferentes do padrão
Principalmente em colunas e barras, faça o eixo começar de ZERO
loadfonts(device = "win")
by_country <- organdata %>% group_by(consent_law, country) %>%
summarize_if(is.numeric, list(mean = mean, sd = sd), na.rm = TRUE) %>%
ungroup()
p <- ggplot(data = by_country,
mapping = aes(x = gdp_mean, y = health_mean))
p + geom_point( color = "coral4") +
geom_text_repel(data = subset(by_country, gdp_mean > 25000),
mapping = aes(label = country),
size = 3,
color = "coral4"
#family="Bookman Old Style"
) +
labs(y = "Gastos com saúde per cap.", x = "PIB per cap." ) +
scale_x_continuous(
labels = number_format(decimal.mark = ",", big.mark = "."),
limits = c(0,NA)
) +
scale_y_continuous(
labels = number_format(decimal.mark = ",", big.mark = "."),
limits = c(0,NA)
) +
theme_minimal(
) +
ggtitle(
label = "Gastos em saúde x PIB"
) +
theme(
#text = element_text(family="Bookman Old Style", color = "coral4"),
text = element_text( color = "coral4"),
axis.text = element_text(colour = "coral4"),
panel.background = element_rect(fill = "beige"),
panel.grid.minor = element_line(colour = "bisque1"),
panel.grid.major = element_line(colour = "bisque3")
)Use small multiples: divida o gráfico em gráficos iguais menores cada um com parte dos dados, seguindo uma regra categórica
#("azure2", "chocolate", "steelblue4", "darkslategray", "slategray3", "slategray4")
loadfonts(device = "win")
p <- ggplot(data = gapminder, mapping = aes(x = year, y = gdpPercap))
p + geom_line(color="chocolate", aes(group = country)) +
geom_smooth(size = 1.1, method = "loess", se = FALSE, color = "tan4") +
scale_y_log10(labels=scales::dollar) +
facet_wrap(~ continent, ncol = 5) +
labs(x = "Year",
y = "GDP per capita",
title = "GDP per capita on Five Continents")+
theme_minimal() +
theme(
text=element_text(family="CMU Serif", color = "darkslategray", size = 8),
axis.text = element_text(colour = "darkslategray"),
panel.background = element_rect(fill = "azure2"),
panel.grid.minor = element_line(colour = "slategray2"),
panel.grid.major = element_line(colour = "slategray2"),
strip.text = element_text(family="CMU Serif", colour = "darkslategray"),
aspect.ratio = 1
)## `geom_smooth()` using formula 'y ~ x'
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
ggplot2 é a biblioteca de visualização de dados moderna mais usada no R.
Muitos dos problemas listados anteriormente são tratados por ela. É até difícil causar alguns dos problemas anteriores, por exemplo as distorções. É preciso muito malabarismo para produzir um gráfico com 2 eixos.
Por outro lado, a biblioteca facilita muito a criação de small multiples, a criação de gráficos personalizados e complexos
GGPLOT é baseada na filosofia de Tufte (Tufte 1973) e Wilkinson (Wilkinson 1999), em que os gráficos são construídos a partir de dois princípios:
Um gráfico é uma construção em camadas
Os dados ganham sentido a partir de mapeamentos a escalas, que são mapeadas a geometrias
A GGPLOT2 parece mais complicada de usar do que funções como a plot() da biblioteca base do R.
Talvez porque as pessoas não são sempre apresentadas às camadas da gramática “Grammar of Graphics” que fundamenta a biblioteca.
As camadas da ggplot
Fonte: (Scavetta 2018)
As camadas da ggplot são as seguintes:
Dados (data)
Aesthetics (aes()): mapeamento das colunas dos dados a escalas
Geometries geoms_(): elementos visuais que mostram os dados nas escalas
Facets: small multiples
Statistics: representações do dados transformadas por operações matemáticas
Coordinates
Themes: “non-data ink,” ou seja, os elementos que não são diretamente plotados pela existência de dados
Fonte: (Scavetta 2018)
Usando as camadas, vamos montando o gráfico
Vamos usar para este exemplo um dataset clássico, de espécimes de flores: iris
## Rows: 150
## Columns: 5
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4,...
## $ Sepal.Width <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7,...
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5,...
## $ Petal.Width <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2,...
## $ Species <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa...
Veja que, pela variável ser representada de forma discreta (uma casa decimal), há sobreposição de pontos. para dar a real impressão de quantos pontos existem, o ideal é inserir um ruído com geom_jitter().
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6) +
facet_grid(. ~ Species)ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6) +
facet_grid(. ~ Species) +
stat_smooth(method = "lm", se = F, col = "red")## `geom_smooth()` using formula 'y ~ x'
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6) +
facet_grid(. ~ Species) +
stat_smooth(method = "lm", se = F, col = "red") +
coord_flip()## `geom_smooth()` using formula 'y ~ x'
wes <- wes_palette("Royal2", 5, "discrete")
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6, color = wes[1]) +
facet_grid(. ~ Species) +
stat_smooth(method = "lm", se = F, col = wes[1]) +
coord_flip() +
theme_minimal() +
theme(
text=element_text( color = wes[1], size = 16),
axis.text = element_text(colour = wes[1]),
panel.background = element_rect(fill = "white"),
panel.grid.minor = element_line(colour = wes[2]),
panel.grid.major = element_line(colour = wes[3]),
panel.border = element_rect(colour = wes[4], fill = NA)
)## `geom_smooth()` using formula 'y ~ x'
gols <- read_csv("dados/football_events/ginf.csv") %>%
select(fthg, ftag) %>%
pivot_longer(cols = everything(), names_to = "casa_fora", values_to = "gols")
ggplot(gols) +
geom_histogram(
aes(x = gols),
binwidth = 1
) +
scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
theme_minimal()result <- fitdistr( gols$gols , densfun = "Poisson" )
pois <- rpois(length(gols$gols ), lambda = result$estimate ) %>%
enframe(value = "gols") %>%
mutate(tipo = "poisson")
real_simulado <- gols %>%
mutate(tipo = "real" ) %>%
bind_rows(pois)
ggplot() +
geom_histogram(
data = real_simulado,
aes(x = gols, group = tipo, fill = tipo, color = tipo),
binwidth = 1,
show.legend = TRUE,
position = "identity",
alpha = 0.3,
size = 2
) +
scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
geom_vline(aes(xintercept = result$estimate)) +
geom_text_repel(
aes(x = result$estimate, y = 8000 ),
nudge_x = 0.45,
label = number(result$estimate, accuracy = 0.01, decimal.mark = "," )
) +
scale_y_continuous(
labels = number_format(big.mark = ".", decimal.mark = ",")
) +
theme_minimal() +
theme(
legend.position = "top"
) +
labs(
fill = "Distribuição",
color = "Distribuição",
x = "\nGols",
y = "Ocorrências\n"
)O gráfico mostrando gráfco com a função de probablidade acumulada empírica mostra propriedades que o histograma não mostra
ggplot(gols, aes(x = gols)) +
stat_ecdf(
geom = "step",
fill = "darkblue",
color = "darkblue",
size = 2
) +
scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
theme_minimal() +
labs(
y = "Quantil\n",
x = "\nGols"
)## Warning: Ignoring unknown parameters: fill
Para variáveis contínuas, faz sentido usar o gráfico de densidade de probabilidade
Nesta transformação, repare o uso da pivot_longer(). Neste caso, o nome das colunas que estão sendo desmobilizadas e transformadas em valores contêm duas informações, separadas por “_”:
values_to)Além disso, uma transformação é aplicada nos nomes das colunas de forma a possibilitar o tratamento no mutate posterior.
ufc_pesos <- read_csv("dados/ufc/data.csv") %>%
filter(weight_class == "Heavyweight") %>%
pivot_longer(
cols = matches("[BR]_.*"),
names_pattern = "(.+?)_(.*)",
names_to = c("lutador", ".value"),
names_transform = list( lutador = ~if_else(.x == "R", "Red", "Blue") )
) %>%
janitor::clean_names() %>%
mutate(
resultado = if_else(lutador == winner, "Vencedor", "Perdedor")
) ggplot(ufc_pesos) +
geom_density(
aes(x = height_cms),
color = "darkblue",
fill = "darkblue",
alpha = 0.4
) +
theme_minimal() +
labs(
y = "",
x = "Altura (cm)"
) +
theme(
axis.text.y = element_blank()
)Mais uma vez a distribuição empírica nos deixa ver mais coisa: os quantis
ggplot(ufc_pesos) +
stat_ecdf(
aes(x = height_cms),
geom = "density" ,
fill = "darkblue",
color = "darkblue",
size = 2,
alpha = 0.5
) +
theme_minimal() +
labs(
x = "\nAltura (cm)",
y = "Quantil\n"
)## Warning: Removed 1 rows containing non-finite values (stat_ecdf).
ggplot(ufc_pesos) +
geom_density(
aes(
x = height_cms,
color = resultado
),
size = 1.5
) +
theme_minimal() +
theme(
legend.position = "top",
axis.text.y = element_blank()
) +
labs(
x = "Altura (cm)",
color = "Resultado do lutador",
fill = "Resultado do lutador",
y = ""
)## Warning: Removed 1 rows containing non-finite values (stat_density).
ggplot(ufc_pesos) +
stat_ecdf(aes(x = height_cms, color = resultado), geom = "density", size = 1.5 ) +
theme_minimal() +
theme(
legend.position = "top"
) +
labs(
x = "Altura (cm)",
color = "Resultado do lutador",
fill = "Resultado do lutador",
y = "Quantil\n"
)## Warning: Removed 1 rows containing non-finite values (stat_ecdf).
O tipo de gráfico gerado por geom_desity_ridges() ajuda a comparar várias distribuições.
ufc_pesos <- read_csv("dados/ufc/data.csv") %>%
select(
weight_class,
R_Height_cms,
B_Height_cms,
Winner
) %>%
transmute(
weight_class,
altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms ),
altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms )
) %>%
pivot_longer(cols = -weight_class, names_to = "resultado", values_to = "altura")
ggplot(ufc_pesos) +
geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
theme_ridges() +
scale_x_continuous(breaks = seq(150,by = 5, to = 220))Legal a ideia, mas não queremos comparar mulheres também
A biblioteca stringr facilita muito lidar com strings. Por exemplo detectar um pedaço de string em outra.
ufc_pesos <- read_csv("dados/ufc/data.csv") %>%
filter(!str_detect(weight_class, "Women")) %>%
select(
weight_class,
R_Height_cms,
B_Height_cms,
Winner
) %>%
transmute(
weight_class,
altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms ),
altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms )
) %>%
pivot_longer(cols = -weight_class, names_to = "resultado", values_to = "altura")
ggplot(ufc_pesos) +
geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
theme_ridges() +
scale_x_continuous(breaks = seq(150,by = 5, to = 220))Melhorou, mas ainda á estranho. Melhor seria se as classes de peso estivessem em ordem
No R, as variáveis categóricas são chamadas de fator. Cada valor possível de um fator é representadas internamente com um identificador numérico, e não com a própria string.
Cada valor possível da variável categórica também é chamada de level
Podemos manipular essas representações facilmente com a forcats.
No caso, queremos ordenar nosso fator com a ordem as classes de peso no UFC
A função fct_relevel possibilita a ordenação manual de uma variável factor
ufc_pesos_ordem <- ufc_pesos %>%
mutate(
weight_class = fct_relevel(weight_class,
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight",
"Catch Weight",
"Open Weight"
)
)
)
ggplot(ufc_pesos_ordem) +
geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
theme_ridges() +
scale_x_continuous(breaks = seq(150,by = 5, to = 220))Os gráficos do tipo scatter plot possibilitam a visualização de relações entre as variáveis
Vamos brincar um pouco com os dados do Banco Mundial
taxa_homicidios <- wb("VC.IHR.PSRC.P5", country = "all", mrv = 10, POSIXct = TRUE ) %>%
select(iso3c, value, date) %>%
group_by(iso3c) %>%
top_n(1, date) %>%
rename(homicidios = value)
gini <- wb("SI.POV.GINI", country = "all", mrv = 10, POSIXct = TRUE ) %>%
select(iso3c, value, date) %>%
group_by(iso3c) %>%
top_n(1, date) %>%
rename(desigualdade = value)
pib_per_capita <- wb("NY.GDP.PCAP.PP.KD", country = "all", mrv = 10, POSIXct = TRUE ) %>%
select(iso3c, value, date) %>%
group_by(iso3c) %>%
top_n(1, date) %>%
rename(riqueza = value)
pib_gini_homi <- taxa_homicidios %>%
inner_join(gini, by = c("iso3c")) %>%
inner_join(pib_per_capita, by = c("iso3c")) %>%
left_join(codelist, by = c("iso3c") ) %>%
select(riqueza, desigualdade, homicidios, continent, region) %>%
pivot_longer(cols = c("riqueza", "desigualdade"), names_to = "tipo_indice", values_to = "indice" ) %>%
mutate(tipo_indice = str_to_title(tipo_indice))ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
geom_point(alpha = 0.4) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth() +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) É possível visualizar os eixos em log na base 10
ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
geom_point(alpha = 0.4) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth() +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() É possível também usar as cores para olhar uma terceira característica
ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
geom_point(alpha = 0.4, aes(color = continent )) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth() +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() +
theme(legend.position = "top")A linha de regressão também pode ser estimada para cada grupo. Neste caso vamos usar o método de regressão linear, pois o LOESS original nao funciona nada bem com muito poucos dados (nem a linear, mas…)
ggplot(pib_gini_homi, aes(y = homicidios, x = indice, color = continent)) +
geom_point(alpha = 0.4 ) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth(se = FALSE, method = "lm") +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() analise_desigualdade <- pib_gini_homi %>%
filter(tipo_indice == "Desigualdade")
ggplot(analise_desigualdade, aes(y = homicidios, x = indice, color = continent)) +
geom_point(alpha = 0.4) +
facet_wrap(~tipo_indice, ncol = 1) +
geom_smooth(se = FALSE, method = "lm") +
theme_minimal() +
labs(x = "Índice de Gini", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() +
facet_wrap(~continent) +
theme(legend.position = "top")ggplot(analise_desigualdade, aes(y = homicidios, x = indice, color = continent)) +
geom_text(aes(label = iso3c), size = 2) +
facet_wrap(~tipo_indice, ncol = 1) +
geom_smooth(se = FALSE, method = "lm") +
theme_minimal() +
labs(x = "Índice de Gini", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() +
facet_wrap(~continent) +
theme(legend.position = "top")Podemos avaliar três variáveis contínuas usando um gráfico
Primeiro vamos colocar em wide as duas variáveis que tínhamos deixado long.
pib_gini_homi_wide <- pib_gini_homi %>%
pivot_wider(names_from = tipo_indice, values_from = indice )
head(pib_gini_homi_wide)## # A tibble: 6 x 6
## # Groups: iso3c [6]
## iso3c homicidios continent region Riqueza Desigualdade
## <chr> <dbl> <chr> <chr> <dbl> <dbl>
## 1 ALB 2.29 Europe Europe & Central Asia 13962. 33.2
## 2 DZA 1.36 Africa Middle East & North Africa 11350. 27.6
## 3 AGO 4.85 Africa Sub-Saharan Africa 6654. 51.3
## 4 ARG 5.32 Americas Latin America & Caribbean 22034. 41.4
## 5 ARM 1.69 Asia Europe & Central Asia 13654. 34.4
## 6 AUS 0.892 Oceania East Asia & Pacific 49756. 34.4
A escala de cor usada abaixo, Viridis, oferece a mesma sensação para colorido e preto e branco e é feita para ajudar daltônicos.
O log foi usado aqui para evitar que os valores extremos dificultem a visualização da escala de cores.
ggplot(pib_gini_homi_wide) +
geom_point(aes(color = homicidios, x = Desigualdade , y = Riqueza )) +
scale_color_viridis_c(trans = "log10", direction = -1, option = "inferno") +
theme_minimal()hdi <- read_csv("dados/hdi/atlas.csv") %>%
rename(municipio = `município` ) %>%
filter(ano == 2010)
violencia <- pdf_text("dados/atlasviolencia/8099-tabelamunicipiostodossite.pdf") %>%
str_split(pattern = fixed("\n")) %>%
unlist() %>%
enframe( value = "linha") %>%
mutate(
UF = str_extract(linha, "[A-Z]{2}") ,
municipio = str_extract(linha, "(?<=[A-Z]{2} ).+?(?=[0-9])"),
numeros = str_extract_all(linha, "(?<= )[0-9 ,.]*")
) %>%
rowwise() %>%
mutate(
numeros = str_flatten(numeros)
) %>%
filter(str_length(numeros)>1) %>%
mutate(
numeros = str_trim(numeros),
municipio = str_trim(municipio) %>% str_to_upper()
) %>%
separate(
col = numeros,
into = c("populacao", "homicidios", "ocultos", "taxa_homicidios"),
sep = "[ ]+"
) %>%
mutate(taxa_homicidios = parse_number(taxa_homicidios, locale = locale(decimal_mark = ",")) )
ufs <- municipios_ajeitado %>%
select(
UF,
cod_municipio
) %>%
mutate(
cod_uf = cod_municipio %/% 100000
) %>%
select(-cod_municipio) %>%
distinct() %>%
filter(!is.na(cod_uf))
hdi %<>% inner_join(ufs, by = c("uf" = "cod_uf"), suffix = c("_old", ""))
hdi_violencia <- hdi %>%
inner_join(violencia, by = c("municipio", "UF") ) %>%
mutate(
regiao = case_when(
uf %/% 10 == 1 ~ "Norte",
uf %/% 10 == 2 ~ "Nordeste",
uf %/% 10 == 3 ~ "Sudeste",
uf %/% 10 == 4 ~ "Sul",
uf %/% 10 == 5 ~ "Centro-Oeste",
TRUE ~ NA_character_
)
) %>%
mutate(
populacao =
parse_number(
populacao,
locale = locale(grouping_mark = ".")
)
)populacao <- hdi_violencia %>%
mutate(
faixa_populacao = case_when(
populacao < 20000 ~ "Menos 20 mil",
populacao < 80000 ~ "20 mil a 80 mil",
populacao < 300000 ~ "80 mil e 300 mil",
TRUE ~ "> 300 mil"
)
)
ggplot(populacao) +
geom_jitter(aes(color = taxa_homicidios, x = idhm_r , y = idhm_e ), alpha = 0.4 ) +
scale_color_gradientn(
values = c(0,0.20,0.4,0.6,0.8, 1),
# values = c(
# 1,
# quantile(
# probs = 0.3,
# hdi_violencia$taxa_homicidios %>% log10(),
# na.rm = TRUE
# ),
# quantile(
# probs = 0.9,
# hdi_violencia$taxa_homicidios %>% log10(),
# na.rm = TRUE
# )
# ),
colors =
c("blue", "blue", "gray50", "yellow", "red", "red"),
na.value = "white",
trans = "log10",
) +
facet_grid(faixa_populacao ~ regiao) +
theme_minimal() +
theme(
legend.position = "top"
) +
labs(
color = "Taxa de Homicídio",
x = "\nIDH (Renda)",
y = "IDH (Educação)\n"
)## Warning: Transformation introduced infinite values in discrete y-axis
ufc_data <- read_csv("dados/ufc/data.csv") %>%
select(weight_class, no_of_rounds )
ufc_raw_data <- read_csv2("dados/ufc/raw_total_fight_data.csv") %>%
select(last_round)
ufc_tudo <- bind_cols(ufc_data, ufc_raw_data) %>%
filter(no_of_rounds == 3) %>%
group_by(weight_class) %>%
mutate(n_classe = n()) %>%
group_by(weight_class, last_round) %>%
summarise(n_round = n(), n_classe = mean(n_classe)) %>%
mutate(frac_lutas = n_round/n_classe ) %>%
ungroup() %>%
filter( weight_class %in%
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight"
)) %>%
mutate(
weight_class = fct_relevel(weight_class,
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight"
)
)
)O gráfico gerado por geom_tile() ajuda a visualizar duas variáveis discretas ordinais e uma contínua.
ggplot(ufc_tudo, aes(x = last_round, y = weight_class)) +
geom_tile(aes(fill = frac_lutas )) +
geom_text(aes(label = percent(frac_lutas)) ) +
scale_fill_gradient(low = "white", high = "darkgreen", labels = percent ) +
theme_minimal() hdi_desc <- read_csv("dados/hdi/desc.csv")
hdi_tudo <- read_csv("dados/hdi/atlas.csv") %>%
mutate(
regiao = case_when(
uf %/% 10 == 1 ~ "N",
uf %/% 10 == 2 ~ "NE",
uf %/% 10 == 3 ~ "SE",
uf %/% 10 == 4 ~ "S",
uf %/% 10 == 5 ~ "CO",
TRUE ~ NA_character_
)
) %>%
select(
espvida,
anos_estudo = e_anosestudo,
gini,
pind,
pren10ricos,
prentrab,
t_banagua,
regiao
)
ggpairs(hdi_tudo) +
theme_minimal()A função geom_col() é melhor para poucos períodos
Pelo mapeamentos no eixo x e y, haveria valores no mesmo lugar do eixo cartesiano.
Para sanar esse conflito, usamos o position - parâmetro da geom_col().
Pata o caso em que queremos comparar valores relativamente a cada período, devemos empilhar as categorias de forma que a coluna tenha o mesmo tamanho a cada período.
Para isso, usamos o valor fill para o parâmetro position
jogos <- read_csv("dados/football_events/ginf.csv") %>%
mutate(
resultado = case_when(
fthg > ftag ~ "Casa",
fthg < ftag ~ "Fora",
TRUE ~ "Empate"
)
) %>%
count(season, country, resultado)
ggplot(jogos) +
geom_col(aes(y = n, fill = resultado, x = as.factor(season)), position = "fill") +
facet_wrap(~country) +
theme_minimal() +
scale_fill_manual(values = wes_palette("Royal2")) +
scale_y_continuous(labels = percent) +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top"
# text = element_text(family = "Rockwell Condensed")
) +
labs(x = "Temporada", y = "% Vitórias", fill = "Resultado")Para compararmos de forma que o total de cada período apareça, em valores absolutos, devemos usar o valor stack
jogos <- read_csv("dados/football_events/ginf.csv") %>%
filter(season != 2017) %>%
select(season, country, Fora = ftag, Casa = fthg) %>%
pivot_longer(cols = c(Fora, Casa), names_to = "Time", values_to = "Gols" ) %>%
group_by(season, country, Time) %>%
summarise(Gols = sum(Gols))
ggplot(jogos) +
geom_col(aes(y = Gols, fill = Time, x = as.factor(season)), position = "stack") +
facet_wrap(~country) +
theme_minimal() +
scale_fill_manual(values = wes_palette("Royal2")) +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top"
# text = element_text(family = "Rockwell Condensed")
) +
labs(x = "Temporada", y = "Gols", fill = "Time")Para muito períodos o ideal é usar áreas empilhadas, com a função geom_area()
Repare que há muitos países, ou seja, muitas categorias.
Usando a biblioteca forcats e sua função fct_lump() conseguimos criar uma categoria de “Outros”
gdps <- wb(indicator = "NY.GDP.MKTP.PP.KD", country = "countries_only")
gdps_tratado <- gdps %>%
mutate(date = as.integer(date) ) %>%
group_by(country) %>%
mutate(ultimo_gdp = last(value, order_by = date)) %>%
ungroup() %>%
mutate(country = fct_lump(country, n = 7, w = ultimo_gdp, other_level = "Outros")) %>%
group_by(country, date) %>%
summarise(value = sum(value),
ultimo_gdp = sum(ultimo_gdp)) %>%
ungroup() %>%
mutate(country = fct_reorder(country, ultimo_gdp ))
ggplot(gdps_tratado ) +
geom_area(aes(x = date, y = value, fill = country), position = "fill") +
theme_minimal() +
scale_fill_brewer(palette = "Accent") +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top",
# text = element_text(family = "CMU Classical Serif")
) +
labs(x = "Ano", y = "PIB PPP", fill = "País")gdps_tratado <- gdps %>%
mutate(date = as.integer(date) ) %>%
group_by(country) %>%
mutate(ultimo_gdp = last(value, order_by = date)) %>%
ungroup() %>%
mutate(country = fct_lump(country, n = 7, w = ultimo_gdp, other_level = "Outros")) %>%
group_by(country, date) %>%
summarise(value = sum(value),
ultimo_gdp = sum(ultimo_gdp)) %>%
ungroup() %>%
mutate(country = fct_reorder(country, ultimo_gdp ))## `summarise()` regrouping output by 'country' (override with `.groups` argument)
ggplot(gdps_tratado ) +
geom_area(aes(x = date, y = value, fill = country) ) +
theme_minimal() +
scale_fill_brewer(palette = "Accent") +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top",
# text = element_text(family = "CMU Serif Extra")
) +
labs(x = "Ano", y = "PIB PPP", fill = "País")Aqui cabe uma torta se houver poucas categorias
Para gerar um gráfico de torta é preciso gerar um gráfico de barras e tornar a coordenada polar com uso de coord_polar
gdps_tratado <- gdps %>%
mutate(date = as.integer(date)) %>%
filter(date == max(date) ) %>%
mutate(country = fct_lump(country, n = 1, w = value, other_level = "Resto")) %>%
group_by(country) %>%
summarise(PIB = sum(value))## `summarise()` ungrouping output (override with `.groups` argument)
ggplot(gdps_tratado, aes(x = "", y = PIB, fill = country)) +
geom_bar( width = 1, stat = "identity", color = "white") +
coord_polar("y", start = 0) +
scale_fill_brewer(palette = "Accent") +
theme_void()gdps_tratado <- gdps %>%
mutate(date = as.integer(date)) %>%
filter(date == max(date) ) %>%
mutate(country = fct_lump(country, n = 50, w = value, other_level = "Resto")) %>%
group_by(country) %>%
summarise(PIB = sum(value))## `summarise()` ungrouping output (override with `.groups` argument)
ggplot(gdps_tratado, aes( fill = country, area = PIB, label = country)) +
geom_treemap() +
geom_treemap_text(grow = TRUE) +
theme(
legend.position = "none"
)Gráficos animados devem ser usados com moderação.
A gganimate ajuda a faxer animações depois que você tiver pensado duas vezes e sucumbido a elas.
Ela é uma camada adicionada a um gráfico ggplot que define um particionamento do tibble referenciado pelo gráfico de forma a gerar múltiplos frames.
A função transition_time() define o campo que particiona o gráfico.
library(gapminder)
plot <- ggplot(gapminder, aes(gdpPercap, lifeExp, size = pop, colour = country)) +
geom_point(alpha = 0.7, show.legend = FALSE) +
scale_colour_manual(values = country_colors) +
scale_size(range = c(2, 12)) +
scale_x_log10() +
facet_wrap(~continent) +
# Here comes the gganimate specific bits
labs(title = 'Year: {frame_time}', x = 'GDP per capita', y = 'life expectancy') +
transition_time(year) +
ease_aes('linear')
animate(plot, nframes = 500, fps = 15, width = 400, height = 400, renderer = gifski_renderer())
anim_save(filename = "imagens/animate_1.gif")Vocês provavelmente já se depararam com a famigerada “)bar racer a_r(.”
Veja como construir sua própria bar race
fundos_escolhidos <- tribble(
~cnpj, ~nome,
"07.455.507/0001-89", "Verde",
"12.798.221/0001-36", "Nimitz",
"17.087.932/0001-16", "Maraú",
"31.666.755/0001-53", "Legacy",
"23.884.632/0001-60", "Adam",
"11.419.627/0001-06", "Itaú"
)
cores <- c(
"Nimitz" = "seashell4",
"Verde" = "seagreen",
"Maraú" = "darkred",
"Itaú" = "darkorange3",
"Legacy" = "lightgoldenrod4",
"Adam" = "lightblue"
)
dados_fundos_escolhidos <- todos_os_fundos %>%
inner_join(
fundos_escolhidos,
by = c("CNPJ_FUNDO" = "cnpj")
) %>%
filter(
!is.na(VL_PATRIM_LIQ) & VL_PATRIM_LIQ > 1000000000
)
dados_fundos_escolhidos_rank <- dados_fundos_escolhidos %>%
group_by(DT_COMPTC) %>%
mutate(
casas = log10(max(VL_PATRIM_LIQ)) %>% floor(),
VL_PATRIM_LIQ = VL_PATRIM_LIQ / 10^(casas - 9),
rank = rank(VL_PATRIM_LIQ),
rank = (6 - max(rank)) + rank,
pl = VL_PATRIM_LIQ/max(VL_PATRIM_LIQ)
) %>%
select(
casas,
everything()
) %>%
ungroup() %>%
mutate(
PL_texto = number(VL_PATRIM_LIQ / 1000000, accuracy = 1, big.mark = ".", decimal.mark = ","),
label = str_glue("{nome}: {PL_texto}")
)
my_plot <- ggplot(dados_fundos_escolhidos_rank, aes(group = nome, y = nome, x = rank, label = nome, fill = nome) ) +
geom_tile(
aes(
y = pl/2,
height = pl,
width = 0.9,
fill = nome
)
) +
geom_text(
aes(
y = pl,
label = label,
hjust = ifelse(pl > 0.5, 1, 0),
size = 3
)
) +
coord_flip(clip = "off", expand = FALSE) +
scale_fill_manual(values = cores) +
# Here comes the gganimate specific bits
labs(title = 'Data: {frame_time}', x = '', y = '') +
transition_time(DT_COMPTC) +
ease_aes("cubic-in-out") +
theme_minimal() +
theme(
plot.title = element_text(color = "#01807e", face = "bold", hjust = 0, size = 30),
axis.ticks.y = element_blank(), #removes axis ticks
axis.text.y = element_blank(), #removes axis ticks
axis.ticks.x = element_blank(), #removes axis ticks
axis.text.x = element_blank(), #removes axis ticks
panel.grid.major = element_blank(), #removes grid lines
panel.grid.minor = element_blank(), #removes grid lines
legend.position = "none",
)
animate(my_plot, nframes = 500, fps = 15, width = 400, height = 400, renderer = gifski_renderer())
anim_save(filename = "imagens/bar_race.gif")A biblioteca sf tem várias funções para lidar com dados geoespaciais.
A função st_read lê vários formatos de dados geoespaciais para um objeto da classe sf, ou simple features.
sf é uma classe que extende a classe data.frame contendo uma coluna de dados geoespaciais.
Cada linha do data.frame é relacionada a um dado geoespacial, que pode ser um ponto, um polígono etc.
## sf [180,674 x 20] (S3: sf/tbl_df/tbl/data.frame)
## $ X : num [1:180674] -39.9 -40.6 -37.1 -48.4 -48.4 ...
## $ Y : num [1:180674] -18.72 -19.38 -11.04 -1.36 -1.38 ...
## $ USER_IdeAg: num [1:180674] 380 381 6587 371 371 ...
## $ USER_SigAg: chr [1:180674] "EDP ES" "ELFSM" "ESE" "CELPA" ...
## $ USER_CodUn: chr [1:180674] "196457" "33243" "703489" "7744250" ...
## $ USER_CodGD: chr [1:180674] "GD.ES.000.034.069" "GD.ES.000.094.399" "GD.SE.000.093.608" "GD.PA.000.108.959" ...
## $ USER_DscCl: chr [1:180674] "Residencial" "Rural" "Residencial" "Residencial" ...
## $ USER_DscGr: chr [1:180674] "B1" "B2" "B1" "B1" ...
## $ USER_DscMo: chr [1:180674] "Autoconsumo remoto" "Autoconsumo remoto" "Geracao na propria UC" "Geracao na propria UC" ...
## $ USER_QtdUC: num [1:180674] 2 2 1 1 1 1 1 1 1 1 ...
## $ USER_CodMu: num [1:180674] 3204906 3203353 2800308 1500800 1500800 ...
## $ USER_NomMu: chr [1:180674] "São Mateus" "Marilândia" "Aracaju" "Ananindeua" ...
## $ USER_SigUF: chr [1:180674] "ES" "ES" "SE" "PA" ...
## $ USER_NomRe: chr [1:180674] "Sudeste" "Sudeste" "Nordeste" "Norte" ...
## $ USER_SigTi: chr [1:180674] "UFV" "UFV" "UFV" "UFV" ...
## $ USER_DscCo: chr [1:180674] "Radiação solar" "Radiação solar" "Radiação solar" "Radiação solar" ...
## $ USER_MdaPo: num [1:180674] 4.6 3 3.96 3 2.4 ...
## $ USER_NomAg: chr [1:180674] "ESPÍRITO SANTO DISTRIBUIÇÃO DE ENERGIA S.A." "EMPRESA LUZ E FORÇA SANTA MARIA S/A" "ENERGISA SERGIPE - DISTRIBUIDORA DE ENERGIA S.A" "CENTRAIS ELÉTRICAS DO PARÁ S.A. - CELPA" ...
## $ DthConexao: Date[1:180674], format: "2018-07-24" "2019-05-21" ...
## $ geometry :sfc_POINT of length 180674; first list element: 'XY' num [1:2] -39.9 -18.7
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:19] "X" "Y" "USER_IdeAg" "USER_SigAg" ...
enderecos <- tribble(
~nome, ~endereco,
"RB1", "Avenida Rio Branco 1, Centro, Rio de Janeiro",
"Marques dos Reis", "Praça Pio X 54, Centro, Rio de Janeiro"
)
enderecos_com_geo <- enderecos %>%
mutate_geocode(location = endereco)## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Avenida+Rio+Branco+1,+Centro,+Rio+de+Janeiro&key=xxx
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Pra%C3%A7a+Pio+X+54,+Centro,+Rio+de+Janeiro&key=xxx
enderecos_id <- enderecos_com_geo %>%
mutate(id = row_number())
comb_1 <- enderecos_id %>%
rename_all(function(x){str_glue("{x}_1")})
comb_2 <- enderecos_id %>%
rename_all(function(x){str_glue("{x}_2")})
ponto_4326 <- function (lat, lon){
st_point(c(lat, lon)) %>%
st_sfc(crs = 4326)
}
enderecos_combinados <- comb_1 %>%
crossing(comb_2) %>%
filter(id_1 < id_2) %>%
mutate(
ponto_1 = map2(.x = lon_1, .y = lat_1, .f = ponto_4326 ) ,
ponto_2 = map2(.x = lon_2, .y = lat_2, .f = ponto_4326 ) ,
dist = map2(.x = ponto_1, .y = ponto_2, st_distance )
)A biblioteca ggmap deixa que usemos um background em um plot de mapa.
Ele funciona em camadas como o ggplot, e aceita as mesmas funcionalidades para plotar os dados, usando o eixo x como latitude e o eixo y como longitude
No código abaixo, usamos get_map para pegar o fundo, passando os parâmetros desejados. Existem várias fontes de dados e vários tipos de fundo
O tema incluído por theme_inset é mais adequado a mapas: não mostra os eixos, com ticks e rótulos, por exemplo.
As coordenadas são mantidas fixas, também, com uso de coord_fixed
map <- get_map(location = "Brazil", zoom = 4, source = "google", maptype = "toner-lite")
ggmap(map) +
coord_fixed() +
theme_inset()Podemos plotar os objetos como pontos no mapa, por exemplo, usando geom_point()
ggmap(map) +
geom_point(
data = gd,
alpha = 0.01,
size = 0.5,
aes(
x = X,
y = Y
),
color = "darkblue"
) +
theme_inset()Usando stat_density_2d, podemos plotar
## maptype = "toner-lite" is only available with source = "stamen".
## resetting to source = "stamen"...
## Source : https://maps.googleapis.com/maps/api/staticmap?center=Brazil&zoom=4&size=640x640&scale=2&maptype=terrain&key=xxx
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Brazil&key=xxx
## Source : http://tile.stamen.com/toner-lite/4/4/7.png
## Source : http://tile.stamen.com/toner-lite/4/5/7.png
## Source : http://tile.stamen.com/toner-lite/4/6/7.png
## Source : http://tile.stamen.com/toner-lite/4/4/8.png
## Source : http://tile.stamen.com/toner-lite/4/5/8.png
## Source : http://tile.stamen.com/toner-lite/4/6/8.png
## Source : http://tile.stamen.com/toner-lite/4/4/9.png
## Source : http://tile.stamen.com/toner-lite/4/5/9.png
## Source : http://tile.stamen.com/toner-lite/4/6/9.png
ggmap(ggmap = map) +
stat_density_2d(
aes(x = X, y = Y, fill = ..level.. ),
bins = 70,
geom = "polygon",
data = gd ,
alpha = .2
) +
scale_fill_gradientn(colors = brewer.pal(10, "YlOrRd")) +
theme_inset()## Warning in brewer.pal(10, "YlOrRd"): n too large, allowed maximum for palette YlOrRd is 9
## Returning the palette you asked for with that many colors
Aqui apresentaremos RUDIMENTOS de Aprendizado Estatístico com o objetivo de desmistificar alguns conceitos e apresentar algumas considerações que nos ajudarão a iniciar um processo deste tipo.
Pesquisando imagens no Google vemos as percepções a respeito de “Machine Learning” e “Statistical Learning”
(Mello 2019)
Temos acesso a um pedaço dos dados existentes na realidade.
É o material que temos para nos ajudar a especificar e treinar o nosso modelo.
Podemos adotar como premissa que existe uma função que traduz como o universo funciona. Esta função determina uma estimativa para o valor \(y\) da variável dependente (ou variável de saída) dado que variáveis independentes (ou variáveis de entrada) assumem valores \(x\).
\(x\) é um vetor de tamanho \(p\). \(p\) é o número de variáveis dependentes que afetam a saída. Esta função não determina perfeitamente \(y\), então temos sempre um \(\epsilon\), um ruído que pode ser devido a efeitos que não podem ser mensurados ou não estão disponíveis para uso.
\[y = f(x) + \epsilon\]
Onde:
\(f(x)\) é uma função desconhecida e \(E(\epsilon) = 0\),
\(\epsilon\) é independente de \(x\) e \(Var(\epsilon) = \sigma^2\) é constante em relação a \(x\).
Nós não temos acesso a esta \(f()\) que representa a VERDADE, mas tentamos estimar uma \(\hat{f}()\).
Nosso modelo de Aprendizado Estatístico é uma fábrica de \(\hat{f}()\) que tenta aproximar a \(f()\) verdadeira.
Existem fábricas para todos os gostos. Desde fábricas que produzem funções simples e inteligíveis, que dependem de poucos parâmetros até funções sem uma forma bem conhecida.
Esta fábrica de \(\hat{f}()\) recebe como insumo variáveis retiradas do ambiente. Essas variáveis podem ser tratadas de forma a deixar as coisas mais fáceis para o modelo, de modo a deixá-lo na cara do gol como Gérson fazia na copa de 70. Esta etapa de preparar as entradas se chama Feature Engineering. Para fazer isso, usamos tudo que aprendemos sobre manipulação de dados.
Modelos paramétricos têm uma forma funcional pré-definida. Partindo dessa forma funcional, escolhe-se um valor para os parâmetros de forma a minimizar alguma função de penalização dentro do conjunto que está sendo usado para treinamento.
Um exemplo é a regressão linear.
\[income = \beta_0 + \beta_{education} education + \beta_{seniority} seniority \]
Neste exemplo assumimos que a relação entre a entrada e os dados de saída é linear, ou seja, observa-se um aumento constante na variável de saída a cada unidade de aumento da variável de entrada ao longo de todo todo o intervalo da variável de entrada.
Já os modelos não-paramétricos não consideram nenhuma premissa a respeito da função \(f()\) a ser estimada.
Os modelos não-paramétricos se adaptam muito melhor aos dados. Mas às vezes se adaptam bem DEMAIS…
Com o(s) modelo(s) treinado(s), faremos predições a respeito de outros dados futuros ou cuja variável dependente ainda desconhecemos.
Podemos também fazer inferências a respeito de como as variáveis independentes afetam a variáveis dependente.
A predição é o ato de descobrir \(y\) a partir de \(x\).
A inferência é o entendimento de como valores diferentes de x afetam y
Podemos dividir o erro de previsão em redutível e irredutível.
Temos que:
\[\hat{y} = \hat{f}(x) \]
Então:
\[E(y - \hat{y})^2 = [f(x) - \hat{f}(x)]^2 + Var(\epsilon)\]
Onde
\([f(x) - \hat{f}(x)]^2\) é um erro redutível: podemos sempre estimar uma \(\hat{f}()\) melhor
e
\(Var(\epsilon)\) é um erro irredutível: não importa quão bem estimemos \(\hat{f}()\) o resíduo \(\epsilon\) está lá
(Hastie, Tibshirani, and Friedman 2001)
Podemos dividir o erro redutível em viés e variância.
Se rodarmos \(\hat{f}(x_0)\) treinando o modelo com vários conjuntos de treinamento e testássemos ele com uma entrada específica qualquer para a qual \(y_0 = f(x_0) + \epsilon\):
\[E[y_o - \hat{f}(x_0)] = Var(\hat{f}(x_0)) + [Bias(\hat{f}(x_0))]^2 + Var(\epsilon)\]
\(Var(\hat{f}(x_0)\) é a variância do modelo: o quanto ele muda quando treinamos ele com outro conjunto de treinamento
\([Bias(\hat{f}(x_0))]^2\) é o erro introduzido por tentarmos aproximar um problema complexo da vida real com uso de um modelo simplificado, também chamado de viés do modelo.
Os modelos variam muito com relação a “transparência da caixa” e essa relação tem correlação com viés e variância.
Existe um trade-off entre viés e variância, que são medidas que a gente não observa diretamente.
Há uma relação inversamente proporcional (grosso modo) entre:
Flexibilidade do modelo, ou poder de se adaptar a uma relação complexa entre entrada e saída. Modelos com maior flexibilidade tendem a ter viés menor e variância maior.
Interpretabilidade, ou facilidade de o usuário do modelo (piloto) inferir as relações entre as variáveis independentes e as variáveis dependentes. Modelos mais fáceis de interpretar tendem a ser menos poderosos portanto tem viés maior. Mas também variam menos.
(James et al. 2013)
Além da falta de interpretabilidade, outra questão pode contar contra os modelos mais complexos, com mais poder de aprender nuances sobre os dados: eles podem aprender características específicas da amostra de treinamento que não podem ser generalizados para o resto dos dados.
Existe um ponto ótimo no poder de aprendizado para que o modelo atinja a melhor performance possível na população como um todo.
Além desse ponto ótimo, onde o erro no conjunto de teste é mínimo, o erro no conjunto de treinamento continua a diminuir, mas o aumento do erro no conjunto de teste mostra que ele perde o poder de generalização. Neste região diz-se que está acontecendo o fenômeno do overfitting.
(Hastie, Tibshirani, and Friedman 2001)
O procedimento mais usado para preparar um modelo para a vida real envolve dividir a nossa base em três pedaços:
Dados de treinamento: que serão usados para treinar o modelo e não servem para avaliar seu desempenho
Dados de validação: servem para comparar modelos de forma que nos possibilite escolher UM modelo com UMA especificação de hiperparâmetros (número de nós da rede neural, sensibilidade da penalização pelo tamanho e número dos coeficientes nas regressões lasso e ridge etc. )
Dados de teste: são separados logo no início do processo, e só são usados numa avaliação final do ÚNICO modelo escolhido.
Model Selection: estimação da performance de vários modelos a fim de escolher o melhor
Model Assessment: tendo escolhido o melhor modelo, estimação do erro de generalização em novos dados
Uma forma de usar todos os dados (EXCETO OS DADOS DE TESTE) tanto para treinar quando para validar é revezando que partes dos dados estão sendo usados para treinamento e para validação.
Este esquema é chamado de k-fold cross validation
Podemos dividir os dados em \(k\) partes. Cada vez uma parte é escolhida para ser a validação.
(Robot 2019)
Primeiro podemos fazer nossa exploração preliminar.
Como exemplo, podemos executar um pequeno tratamento das features usando a função pivot_longer() de um jeito diferente, só possível na última versão da biblioteca tidyr.
ufc_data <- read_csv("dados/ufc/data.csv")
ufc_raw_data <- read_csv2("dados/ufc/raw_total_fight_data.csv")
ufc_raw_fighter <- read_csv2("dados/ufc/raw_fighter_details.csv")
ufc_raw_pre <- read_csv2("dados/ufc/preprocessed_data.csv")
ufc_cada_lutador <- bind_cols(ufc_data, ufc_raw_data ) %>%
rename_with(
~str_remove(.x,"[\\.]{3}[0-9]")
) %>%
pivot_longer(
cols = matches("[BR]_.*"),
names_to = c("lutador",".value") ,
names_sep = "_"
) %>%
mutate(
ganhou = str_sub(Winner,1,1) == lutador
) %>%
mutate(aprov = wins/(wins+losses) )
ggplot(ufc_cada_lutador) +
geom_boxplot(aes(x = ganhou, y = aprov )) +
facet_wrap(~weight_class) +
theme_minimal()ufc_aprend <- bind_cols(ufc_data, ufc_raw_data) %>%
rename_with(
~str_remove(.x,"[\\.]{3}[0-9]")
) %>%
filter(Winner != "Draw") %>%
mutate(
R_aprov = R_wins / (R_wins + R_losses),
B_aprov = B_wins / (B_wins + B_losses),
) %>%
select(
Winner,
weight_class,
B_Weight_lbs,
B_age,
R_Weight_lbs,
R_age,
B_age,
R_aprov,
B_aprov
) %>%
mutate(
winner_color = Winner,
zebra_blue = if_else(Winner == "Blue",1,0),
red_mais_velho = R_age - B_age,
idade_red = R_age
) %>%
filter(
!is.na(idade_red) & !is.na(red_mais_velho) & !is.na(zebra_blue)
)Como exemplo, podemos observar a relação entre a diferença de idade dos lutadores e o vencedor da luta.
Normalmente o lutador vermelho é o favorito. mas podemos ver que quando o lutador vermelho é velho e mais velho que o azul, ele costuma perder.
ggplot(ufc_aprend ) +
geom_jitter(aes( y = red_mais_velho, x = R_age, color = winner_color ), alpha = 0.15) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
theme_minimal()Podemos tentar uma regressão linear para avaliar isso:
\[ zebra = \alpha + \beta_{idade} idade\_vermelho + \beta_{dif} dif\_vermelho\_mais\_velho + \epsilon \]
Podemos perceber que os betas são significativos.
##
## Call:
## lm(formula = zebra_blue ~ idade_red + red_mais_velho, data = ufc_aprend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6314 -0.3407 -0.2586 0.5898 0.8888
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.051147 0.061605 0.830 0.406
## idade_red 0.009233 0.002089 4.420 1.01e-05 ***
## red_mais_velho 0.010876 0.001651 6.589 4.91e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4609 on 4865 degrees of freedom
## Multiple R-squared: 0.03404, Adjusted R-squared: 0.03365
## F-statistic: 85.73 on 2 and 4865 DF, p-value: < 2.2e-16
broomA biblioteca broom ajuda a extrair informações de vários tipos de modelo com as mesmas funções.
## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.0511 0.0616 0.830 4.06e- 1
## 2 idade_red 0.00923 0.00209 4.42 1.01e- 5
## 3 red_mais_velho 0.0109 0.00165 6.59 4.91e-11
## # A tibble: 6 x 10
## zebra_blue idade_red red_mais_velho .fitted .se.fit .resid .hat .sigma
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 32 1 0.357 0.00806 -0.357 3.06e-4 0.461
## 2 0 31 -1 0.326 0.00819 -0.326 3.16e-4 0.461
## 3 0 35 -1 0.363 0.0146 -0.363 1.00e-3 0.461
## 4 1 29 3 0.352 0.00839 0.648 3.31e-4 0.461
## 5 1 26 -6 0.226 0.0103 0.774 5.03e-4 0.461
## 6 0 28 -5 0.255 0.00973 -0.255 4.46e-4 0.461
## # ... with 2 more variables: .cooksd <dbl>, .std.resid <dbl>
caretA biblioteca caret também ajuda muito.
A função confusionMatrix monta uma matriz de confusão, que mostra estão acontecendo os falsos positivos, falsos negativos, verdadeiros positivos e verdadeiros negativos.
modelo_aumentado_factor <- modelo_aumentado %>%
mutate(
previsao_zebra_blue = as_factor(round(.fitted)),
zebra_blue = as_factor(zebra_blue)
)
confusionMatrix(modelo_aumentado_factor$previsao_zebra_blue, modelo_aumentado_factor$zebra_blue )## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 3227 1520
## 1 53 68
##
## Accuracy : 0.6769
## 95% CI : (0.6635, 0.69)
## No Information Rate : 0.6738
## P-Value [Acc > NIR] : 0.3293
##
## Kappa : 0.035
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.98384
## Specificity : 0.04282
## Pos Pred Value : 0.67980
## Neg Pred Value : 0.56198
## Prevalence : 0.67379
## Detection Rate : 0.66290
## Detection Prevalence : 0.97514
## Balanced Accuracy : 0.51333
##
## 'Positive' Class : 0
##
purrPreparando um grid das variáveis independentes
red_mais_velho <- seq(
from = min(ufc_aprend$red_mais_velho, na.rm = TRUE),
to = max(ufc_aprend$red_mais_velho, na.rm = TRUE),
length.out = 50
) %>%
enframe(name = "1", value = "red_mais_velho")
idade_red <- seq(
from = min(ufc_aprend$idade_red, na.rm = TRUE),
to = max(ufc_aprend$idade_red, na.rm = TRUE),
length.out = 50
) %>%
enframe(name = "2", value = "idade_red")
weight_classes <- ufc_aprend %>%
select(weight_class) %>%
distinct()
grid_previsao <- red_mais_velho %>%
crossing(idade_red) %>%
crossing(weight_classes)
ggplot(grid_previsao %>% filter(weight_class == "Heavyweight")) +
geom_point(aes(x = idade_red, y = red_mais_velho ), size = 0.01) +
theme_minimal()purrRodando o modelo com o grid
previsoes <- augment(modelo_lm, newdata = grid_previsao )
#idade_red + red_mais_velho
ggplot(ufc_aprend ) +
geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
geom_point(
data = previsoes %>% filter(.fitted > 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightblue",
alpha = 0.1
) +
geom_point(
data = previsoes %>% filter(.fitted < 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightcoral",
alpha = 0.1
) +
theme_minimal()purrrPodemos treinar o modelo para cada classe usando group_by + nest() + map().
Estes comandos criam um tibble com uma coluna de tibbles aninhados já separados por classe de peso. A função map() pode, então rodar o modelo para cada classe.
modelo_ufc_por_classe <- ufc_aprend %>%
group_by(weight_class) %>%
nest_legacy() %>%
mutate(
modelo = map(data, ~lm(zebra_blue ~ idade_red + red_mais_velho, data = .x)),
) %>%
select(-data)
ufc_por_classe <- modelo_ufc_por_classe %>%
mutate(
aumentado = map(modelo, augment)
) %>%
unnest(aumentado)
ufc_por_classe_treinamento <- ufc_por_classe %>%
ungroup() %>%
mutate(
previsao_zebra_blue = as_factor(round(.fitted)),
zebra_blue = as_factor(zebra_blue)
)
confusionMatrix(ufc_por_classe_treinamento$previsao_zebra_blue, ufc_por_classe_treinamento$zebra_blue )## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 3128 1376
## 1 152 212
##
## Accuracy : 0.6861
## 95% CI : (0.6729, 0.6991)
## No Information Rate : 0.6738
## P-Value [Acc > NIR] : 0.03414
##
## Kappa : 0.1088
##
## Mcnemar's Test P-Value : < 2e-16
##
## Sensitivity : 0.9537
## Specificity : 0.1335
## Pos Pred Value : 0.6945
## Neg Pred Value : 0.5824
## Prevalence : 0.6738
## Detection Rate : 0.6426
## Detection Prevalence : 0.9252
## Balanced Accuracy : 0.5436
##
## 'Positive' Class : 0
##
grid_previsao_com_modelo <- grid_previsao %>%
group_by(weight_class) %>%
nest() %>%
inner_join(modelo_ufc_por_classe, by = c("weight_class")) %>%
mutate(previsoes = map2(.x = modelo, .y = data, .f = ~augment(x = .x, newdata = .y) )) %>%
unnest_legacy(previsoes)
ggplot(ufc_aprend ) +
geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
geom_point(
data = grid_previsao_com_modelo %>% filter(.fitted > 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightblue",
alpha = 0.1
) +
geom_point(
data = grid_previsao_com_modelo %>% filter(.fitted < 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightcoral",
alpha = 0.1
) +
theme_minimal()set.seed(2512)
ufc_aprend_classes_selec <- ufc_aprend %>%
filter(weight_class %in%
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight"
)
)
index_treino <- createDataPartition(
ufc_aprend_classes_selec$zebra_blue,
p = 0.75,
list = FALSE,
times = 1
)
ufc_aprend_treino <- ufc_aprend_classes_selec %>%
filter(row_number() %in% index_treino) %>%
select(
idade_red,
red_mais_velho,
zebra_blue,
weight_class
)
ufc_aprend_teste <- ufc_aprend_classes_selec %>%
filter(!(row_number() %in% index_treino)) %>%
select(
idade_red,
red_mais_velho,
zebra_blue,
weight_class
)
modelo_ufc_por_classe_treino <- ufc_aprend_treino %>%
group_by(weight_class) %>%
nest_legacy() %>%
mutate(
modelo = map(data, ~lm(zebra_blue ~ idade_red + red_mais_velho, data = .x)),
) %>%
select(-data)
resultados_teste_lm <- ufc_aprend_teste %>%
group_by(weight_class) %>%
nest_legacy() %>%
inner_join(modelo_ufc_por_classe_treino, by = c("weight_class") ) %>%
mutate(
previsao = map2(.x = modelo, .y = data, .f = ~predict.lm(.x, .y) )
) %>%
unnest(cols = c("data","previsao")) %>%
ungroup() %>%
mutate(
actual = if_else(zebra_blue == 1, "B", "R"),
predicted = if_else(previsao > 0.5, "B", "R"),
) %>%
mutate(
actual = fct_relevel(actual, "R", "B"),
predicted = fct_relevel(predicted, "R", "B")
) %>%
select(
actual,
predicted
)
confusionMatrix(resultados_teste_lm$predicted, resultados_teste_lm$actual )## Confusion Matrix and Statistics
##
## Reference
## Prediction R B
## R 733 325
## B 31 38
##
## Accuracy : 0.6841
## 95% CI : (0.6561, 0.7112)
## No Information Rate : 0.6779
## P-Value [Acc > NIR] : 0.3405
##
## Kappa : 0.0814
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.9594
## Specificity : 0.1047
## Pos Pred Value : 0.6928
## Neg Pred Value : 0.5507
## Prevalence : 0.6779
## Detection Rate : 0.6504
## Detection Prevalence : 0.9388
## Balanced Accuracy : 0.5321
##
## 'Positive' Class : R
##
#ALGO ERRADO AO RODAR O MARKDOWN. NÃO SEI
# modelo_ufc_por_classe_treino_gam <- ufc_aprend_treino %>%
# # group_by(weight_class) %>%
# nest(data =c(zebra_blue, idade_red, red_mais_velho ) ) %>%
# mutate(
# modelo = map(data, ~gam(formula = zebra_blue ~ s(idade_red) + s(red_mais_velho), data = .x))
# ) %>%
# select(-data)
#
#
# resultados_teste_gam <- ufc_aprend_teste %>%
# group_by(weight_class) %>%
# nest_legacy() %>%
# inner_join(modelo_ufc_por_classe_treino_gam, by = c("weight_class") ) %>%
# mutate(
# previsao = map2(.x = modelo, .y = data, .f = ~predict(.x, .y) )
# ) %>%
# unnest(cols = c("data","previsao")) %>%
# ungroup() %>%
# mutate(
# actual = if_else(zebra_blue == 1, "B", "R"),
# predicted = if_else(previsao > 0.5, "B", "R"),
# ) %>%
# mutate(
# actual = fct_relevel(actual, "R", "B"),
# predicted = fct_relevel(predicted, "R", "B")
# ) %>%
# select(
# actual,
# predicted
# )
#
# write_rds(resultados_teste_gam, "dados/ufc/resultados_teste_gam.rds")
resultados_teste_gam <- read_rds("dados/ufc/resultados_teste_gam.rds")
confusionMatrix(resultados_teste_gam$predicted, resultados_teste_gam$actual )## Confusion Matrix and Statistics
##
## Reference
## Prediction R B
## R 716 310
## B 48 53
##
## Accuracy : 0.6823
## 95% CI : (0.6543, 0.7095)
## No Information Rate : 0.6779
## P-Value [Acc > NIR] : 0.3884
##
## Kappa : 0.1026
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.9372
## Specificity : 0.1460
## Pos Pred Value : 0.6979
## Neg Pred Value : 0.5248
## Prevalence : 0.6779
## Detection Rate : 0.6353
## Detection Prevalence : 0.9104
## Balanced Accuracy : 0.5416
##
## 'Positive' Class : R
##
Plotando as regiões deste modelo mais complexo
É possível observar que este modelo se adapta com muito mais potência às especificidades dos dados. Ou seja, consegue um viés baixo no conjunto de treinamento.
Ele se adapta tão bem ao conjunto que se apresentarmos um conjunto um pouco diferente, a \(\hat{f}()\) será diferente certamente.
O modelo parece se adaptar demais aos dados, gerando overfitting. Será que ele é melhor do que um modelo simples?
Precisamos de um mecanismo para avaliar o poder de generalização do modelo.
#
# grid_previsao_com_modelo_gam <- grid_previsao %>%
# group_by(weight_class) %>%
# nest() %>%
# inner_join(modelo_ufc_por_classe_treino_gam, by = c("weight_class")) %>%
# mutate(previsoes = map2(.x = modelo, .y = data, .f = ~predict( .x, .y) )) %>%
# select(-modelo) %>%
# unnest(cols = c("previsoes","data")) %>%
# ungroup()
#
#
# write_rds(grid_previsao_com_modelo_gam, "dados/ufc/grid_previsao_com_modelo_gam.RDS")
grid_previsao_com_modelo_gam <- read_rds("dados/ufc/grid_previsao_com_modelo_gam.RDS")
ggplot( bind_rows(ufc_aprend_classes_selec) ) +
geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
geom_point(
data = grid_previsao_com_modelo_gam %>% filter(previsoes > 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightblue",
alpha = 0.1
) +
geom_point(
data = grid_previsao_com_modelo_gam %>% filter(previsoes < 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightcoral",
alpha = 0.1
) +
theme_minimal()Vamos mostrar a biblioteca caret, que facilita o processo de aprendizado de vários modelos com várias especificações de hiperparâmetros, executando o processo de cross-validation.
Primeiro vamos ler uma base com informações e diagnósticos de câncer de mama (UCI 2019)
Já vamos separar os dados de teste
options(OutDec = ",")
dados <- read_csv("dados/diagnostico/wdbc.csv") %>%
select(-"ID number") %>%
rename_all( .funs = ~str_replace(.,"fractal-media ","fractal_") ) %>%
rename_all( .funs = ~str_replace(.,"fractal-pior ","fractal_") ) %>%
rename_all( .funs = ~str_replace(.,"fractal-dv ","fractal_") ) %>%
rename("fractal_dimension-media" = "fractal_dimension-meia") %>%
rename_all( .funs = ~str_replace(.,"-","_") ) %>%
rename_all( .funs = ~str_replace(.," ","_") ) %>%
mutate(Diagnosis = as.factor(Diagnosis)) ## Parsed with column specification:
## cols(
## .default = col_double(),
## Diagnosis = col_character()
## )
## See spec(...) for full column specifications.
A função trainControl() estabelece os conjuntos para cross validation.
Os mesmos conjuntos serão usados para todos os modelos
A função train() permite o estabelecimento de uma métrica para a escolha do melhor valor para os hiperparâmetros (no caso, não existem), e execuções de funções de pré-processamento.
A função retorna um objeto com todas as informações relativas ao treinamento e a melhor especificação do modelo treinada com todos os dados passados.
A saída padrão impressa em problemas de classificação mostra o resultado do melhor modelo nos seus dados de validação dele.
model_logistic <- train(
form = Diagnosis ~ . ,
data = treino,
metric = "ROC",
method = "glm",
trControl = controle_cv,
preProcess = c( "center", "scale"),
family=binomial(link='logit')
)
model_logistic## Generalized Linear Model
##
## 398 samples
## 30 predictor
## 2 classes: 'B', 'M'
##
## Pre-processing: centered (30), scaled (30)
## Resampling: Cross-Validated (4 fold, repeated 10 times)
## Summary of sample sizes: 299, 298, 299, 298, 299, 298, ...
## Resampling results:
##
## ROC Sens Spec
## 0,9549433 0,9515084 0,9121429
A curva ROC (Receiver Operating Characteristics) foi inventada na época da Segunda Guerra Mundial para avaliar se os operadores de radar americanos estavam detectando confiavelmente aeronaves japonesas a partir de sinais de radar.
A curva mostra, para vários thresholds, qual a fração de verdadeiros positivos (ou sensibilidade) e a fração de falsos positivos (fall-out, ou \(1 - especificidade\) ).
Uma métrica numérica que traduz o quão bom um modelo de classificação é consiste na área embaixo desta curva (AUC, Area Under the Curve). Note que quanto mais perto de um essa área, menor a taxa de falsos positivos e maior a sensibilidade
ggplot(model_logistic$pred, aes(m = M, d = obs, color = Resample )) +
geom_roc( labels = FALSE ) +
coord_equal() + style_roc() + ggtitle("ROC", subtitle = "Métricas para diversos thresholds" )+ style_roc()Podemos avaliar a ROC de cada treinamento feito.
result_model_logistic <- model_logistic$resample %>%
select(Resample, ROC) %>%
rename(AUC = ROC)
result_model_logistic %>%
mutate_if(is.numeric, percent) %>%
kable( caption = "\\label{tab_auc_reglog}Métricas para cada Fold da regressão logistica")| Resample | AUC |
|---|---|
| Fold1.Rep01 | 94.3750% |
| Fold2.Rep01 | 96.1319% |
| Fold3.Rep01 | 93.0357% |
| Fold4.Rep01 | 95.6044% |
| Fold1.Rep02 | 96.4955% |
| Fold2.Rep02 | 93.7582% |
| Fold3.Rep02 | 95.1209% |
| Fold4.Rep02 | 95.2679% |
| Fold1.Rep03 | 95.3125% |
| Fold2.Rep03 | 96.6964% |
| Fold3.Rep03 | 94.2418% |
| Fold4.Rep03 | 96.6813% |
| Fold1.Rep04 | 94.2188% |
| Fold2.Rep04 | 95.4286% |
| Fold3.Rep04 | 95.8901% |
| Fold4.Rep04 | 96.5179% |
| Fold1.Rep05 | 98.1250% |
| Fold2.Rep05 | 99.2308% |
| Fold3.Rep05 | 91.9780% |
| Fold4.Rep05 | 96.4286% |
| Fold1.Rep06 | 98.4152% |
| Fold2.Rep06 | 90.0220% |
| Fold3.Rep06 | 97.0769% |
| Fold4.Rep06 | 96.8750% |
| Fold1.Rep07 | 95.0893% |
| Fold2.Rep07 | 98.1099% |
| Fold3.Rep07 | 96.8352% |
| Fold4.Rep07 | 88.0357% |
| Fold1.Rep08 | 96.0000% |
| Fold2.Rep08 | 90.6813% |
| Fold3.Rep08 | 96.6071% |
| Fold4.Rep08 | 95.3571% |
| Fold1.Rep09 | 96.5275% |
| Fold2.Rep09 | 95.1562% |
| Fold3.Rep09 | 95.3407% |
| Fold4.Rep09 | 96.0714% |
| Fold1.Rep10 | 95.2232% |
| Fold2.Rep10 | 98.4396% |
| Fold3.Rep10 | 95.7143% |
| Fold4.Rep10 | 97.6562% |
result_model_logistic %>%
summarise(media = mean(AUC), sd = sd(AUC)) %>%
rename("Média AUC" = media, "Desvio-padrão AUC" = sd) %>%
mutate_if(is.numeric, percent) %>%
kable(caption = "\\label{tab_auc_reglog_grup}Média e desvio-padrão da Métrica AUC para regressão logistica")| Média AUC | Desvio-padrão AUC |
|---|---|
| 95% | 2% |
texreg(model_logistic$finalModel, custom.model.names = c("Regressão logística simples"), caption = "Coeficientes estimados na regressão logística", label = "reg:log", fontsize = "footnotesize")##
## \begin{table}
## \begin{center}
## \begin{footnotesize}
## \begin{tabular}{l c}
## \hline
## & Regressão logística simples \\
## \hline
## (Intercept) & $80,96$ \\
## & $(35914,01)$ \\
## radius\_media & $-264,43$ \\
## & $(767230,04)$ \\
## texture\_media & $41,40$ \\
## & $(21122,91)$ \\
## perimeter\_media & $1183,18$ \\
## & $(1099949,43)$ \\
## area\_media & $-1040,37$ \\
## & $(494524,57)$ \\
## smoothness\_media & $-65,13$ \\
## & $(30919,21)$ \\
## compactness\_media & $-364,13$ \\
## & $(79670,45)$ \\
## concavity\_media & $6,25$ \\
## & $(189198,65)$ \\
## concave\_points\_media & $367,79$ \\
## & $(89897,83)$ \\
## symmetry\_media & $26,33$ \\
## & $(26581,88)$ \\
## fractal\_dimension\_media & $80,71$ \\
## & $(36331,35)$ \\
## radius\_dv & $246,29$ \\
## & $(104667,71)$ \\
## texture\_dv & $14,71$ \\
## & $(31824,49)$ \\
## perimeter\_dv & $73,35$ \\
## & $(67727,22)$ \\
## area\_dv & $-468,01$ \\
## & $(232868,82)$ \\
## smoothness\_dv & $-73,73$ \\
## & $(23297,51)$ \\
## compactness\_dv & $100,76$ \\
## & $(46184,65)$ \\
## concavity\_dv & $-24,23$ \\
## & $(67029,28)$ \\
## concave\_points\_dv & $51,13$ \\
## & $(46579,32)$ \\
## symmetry\_dv & $71,34$ \\
## & $(24137,24)$ \\
## fractal\_dimension\_dv & $-236,13$ \\
## & $(68224,41)$ \\
## radius\_pior & $-817,84$ \\
## & $(427324,69)$ \\
## texture\_pior & $0,44$ \\
## & $(18864,24)$ \\
## perimeter\_pior & $-187,19$ \\
## & $(129752,63)$ \\
## area\_pior & $1626,90$ \\
## & $(633836,01)$ \\
## smoothness\_pior & $56,67$ \\
## & $(25600,89)$ \\
## compactness\_pior & $-114,58$ \\
## & $(56734,42)$ \\
## concavity\_pior & $95,60$ \\
## & $(74790,04)$ \\
## concave\_points\_pior & $-22,85$ \\
## & $(41737,35)$ \\
## symmetry\_pior & $-28,93$ \\
## & $(38034,29)$ \\
## fractal\_dimension\_pior & $203,37$ \\
## & $(76158,01)$ \\
## \hline
## AIC & $62,00$ \\
## BIC & $185,58$ \\
## Log Likelihood & $-0,00$ \\
## Deviance & $0,00$ \\
## Num. obs. & $398$ \\
## \hline
## \multicolumn{2}{l}{\tiny{$^{***}p<0,001$; $^{**}p<0,01$; $^{*}p<0,05$}}
## \end{tabular}
## \end{footnotesize}
## \caption{Coeficientes estimados na regressão logística}
## \label{reg:log}
## \end{center}
## \end{table}
Uma outra forma de evitar que muitas variáveis sejam usadas no modelo é aplicar uma penalidade diminuir a variância do modelo. É isso que as regressões do tipo Ridge e Lasso fazem.
O modelo mostrado nessa seção é rodado com várias parametrizações. Este modelo, chamado Elastic Net, conjuga a penalização do tipo Ridge com a penalização do tipo Lasso modificando a função de penalização da regressão, que originalmente é o erro quadrático:
\[RSS = \sum_{i = 1}^{n} ( y_i - \beta_0 - \sum_{j=1}^{p}\beta_j x_{ij})^2 \]
Para a regressão Ridge, os coeficientes são penalizados de forma quadrática. Isso diminui a variância do modelo mas não diminui tanto a diminuição do número de coeficientes diferentes de 0:
\[Loss_{Ridge} = RSS + \lambda \sum_{j=1}^{p}\beta_j^2 \]
Para a regressão Lasso, os coeficientes são penalizados pelo seu valor absoluto. Isso diminui a variância do modelo E diminui o número de coeficientes diferentes de 0, favorecendo a interpretabilidade
\[Loss_{Lasso} = RSS + \lambda \sum_{j=1}^{p} \left| \beta_j \right| \]
Cada conjunto de hiperparâmetros do modelo Elastic Net rodado nesta possui dois valores. Um dos valores, \(\alpha\) define que peso deve ser dado a cada tipo de penalização Ridge ou Lasso, onde \(\alpha = 0\) é Ridge puro e \(\alpha = 1\) é Lasso puro. O outro parâmetro, \(\lambda\), regula a intensidade da penalização.
Os resultados para várias parametrizações é mostrado abaixo.
Veja como podemos passar um grid de hiperparâmetros e avaliar esses resultados.
set.seed(13)
model_net <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "glmnet",
trControl = controle_cv,
preProcess = c( "center", "scale"),
tuneGrid = expand.grid(
alpha = c(0,0.5,0.75,0.1,0.125,0.15,0.25,0.5, 1),
lambda = 0:10/500
)
)
plot(model_net)Os resultados de cada parametrização são mostrados em ordem decrescente de AUC.
resultados_net <- model_net$resample %>%
group_by(alpha, lambda) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media) ## `summarise()` regrouping output by 'alpha' (override with `.groups` argument)
## # A tibble: 20 x 4
## # Groups: alpha [6]
## alpha lambda `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl> <dbl> <dbl>
## 1 1 0.02 0.994 0.00746
## 2 1 0.014 0.994 0.00762
## 3 1 0.018 0.994 0.00756
## 4 1 0.016 0.994 0.00761
## 5 1 0.012 0.994 0.00776
## 6 1 0.01 0.994 0.00782
## 7 1 0.008 0.994 0.00777
## 8 0.75 0.006 0.994 0.00790
## 9 0.1 0.01 0.994 0.00715
## 10 0.25 0.01 0.994 0.00732
## 11 0.125 0.018 0.994 0.00732
## 12 0.15 0.008 0.994 0.00720
## 13 0.1 0.02 0.994 0.00730
## 14 1 0.006 0.994 0.00784
## 15 0.1 0.016 0.994 0.00724
## 16 0.1 0.018 0.994 0.00726
## 17 0.15 0.014 0.994 0.00728
## 18 0.1 0.012 0.994 0.00716
## 19 0.75 0.004 0.994 0.00766
## 20 0.125 0.008 0.994 0.00721
Uma outra forma de reduzir a dimensionalidade do probelam é usar o método PCA (Principal Component Analysis) para transformar as variáveis em componentes principais em menor número mas com boa parte da informação original.
Veja como estabelecemos esse pré-processamento.
set.seed(13)
model_net_pca <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "glmnet",
trControl = controle_cv,
preProcess = c( "center", "scale", "pca"),
tuneGrid = expand.grid(
alpha = c(0,0.5,0.75,0.1,0.125,0.15,0.25,0.5, 1),
lambda = 0:10/500
)
)
plot(model_net_pca)resultados_net_pca <- model_net_pca$resample %>%
group_by(alpha, lambda) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)## `summarise()` regrouping output by 'alpha' (override with `.groups` argument)
## # A tibble: 88 x 4
## # Groups: alpha [8]
## alpha lambda `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl> <dbl> <dbl>
## 1 1 0.01 0.993 0.00664
## 2 1 0.012 0.993 0.00658
## 3 1 0.014 0.993 0.00648
## 4 1 0.008 0.993 0.00692
## 5 0.75 0.016 0.993 0.00709
## 6 0.75 0.012 0.993 0.00734
## 7 0.75 0.014 0.993 0.00719
## 8 0.75 0.006 0.993 0.00754
## 9 0.75 0.008 0.993 0.00758
## 10 1 0.006 0.993 0.00703
## # ... with 78 more rows
A árvore de decisão particiona o espaço formado pelas variáveis explicativas em subespaços baseando-se na “pureza” desses subespaços com relação à variável dependente.
model_tree <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "rpart",
trControl = controle_cv,
preProcess = c( "center", "scale")
)
resultados_tree <- model_tree$resample %>%
group_by(cp) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 3 x 3
## cp `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl> <dbl>
## 1 0.0214 0.942 0.0273
## 2 0.0500 0.932 0.0293
## 3 0.829 0.773 0.193
A forma como a árvore de decisão é criada faz com que ela tenha muita variância.
Cada decisão de particionamento é tomada a partir de características que podem ser muito específicas aos dados de treinamento.
Uma ideia usada nas Random Forests é criar conjuntos de treinamento criados a partir do conjunto original, mas retirando amostras de mesmo tamanho COM REPOSIÇÃO. Além disso, cada vez que a árvore é particionada, a partição só pdoe acontecer em \(m\) das variáveis explicativas.
O resultado final é uma média do resultado dessas várias árvores
Estas duas mudanças fazem com que o modelo tenha uma variância muito menor do que as árvore de decisão simples
model_ranger <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "ranger",
trControl = controle_cv,
preProcess = c( "center", "scale"),
verbose = TRUE,
tuneGrid = expand.grid(
mtry = seq(2, by = 2, to = 20),
splitrule = c("gini", "extratrees"),
min.node.size = 1
)
)
plot(model_ranger)resultados_ranger <- model_ranger$resample %>%
group_by(mtry, splitrule) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)## `summarise()` regrouping output by 'mtry' (override with `.groups` argument)
## # A tibble: 20 x 4
## # Groups: mtry [10]
## mtry splitrule `Média AUC` `Desvio-padrão AUC`
## <dbl> <fct> <dbl> <dbl>
## 1 2 extratrees 0.992 0.00606
## 2 4 extratrees 0.992 0.00647
## 3 8 extratrees 0.992 0.00682
## 4 6 extratrees 0.991 0.00686
## 5 16 extratrees 0.991 0.00778
## 6 10 extratrees 0.991 0.00759
## 7 12 extratrees 0.990 0.00808
## 8 14 extratrees 0.990 0.00847
## 9 20 extratrees 0.990 0.00848
## 10 18 extratrees 0.990 0.00913
## 11 2 gini 0.989 0.00933
## 12 6 gini 0.988 0.0104
## 13 18 gini 0.988 0.0101
## 14 16 gini 0.988 0.0101
## 15 10 gini 0.988 0.0104
## 16 4 gini 0.988 0.0105
## 17 20 gini 0.988 0.0101
## 18 8 gini 0.987 0.0106
## 19 14 gini 0.987 0.0104
## 20 12 gini 0.987 0.0108
Podemos montar um modelo ensemble, onde cada modelo original vota e a escolha final recai sobre a classificação que receber mais votos
pred_logistic <- model_logistic$pred %>%
semi_join(model_logistic$bestTune) %>%
mutate(modelo = "Logística")
pred_net_pca <- model_net_pca$pred %>%
semi_join(model_net_pca$bestTune) %>%
mutate(modelo = "Lasso-Ridge PCA")
pred_ranger <- model_ranger$pred %>%
semi_join(model_ranger$bestTune) %>%
mutate(modelo = "Floresta")
pred_ensemble_todos <- bind_rows(pred_logistic, pred_net_pca, pred_ranger)
pred_ensemble_mediana <- pred_ensemble_todos %>%
separate(Resample, c("Fold", "Repeat")) %>%
group_by(Repeat, rowIndex) %>%
summarise(M = median(M), B= median(B), n = n(), obs = unique(obs))
ret_roc <- function(data)
{
roc(data$obs,data$M)
}
ret_roc_auc <- function(data)
{
unlist(as.numeric(roc(data$obs,data$M)$auc))
}
pred_ensemble_mediana_foldado <- model_ranger$pred %>%
semi_join(model_ranger$bestTune) %>%
separate(Resample, c("Fold", "Repeat")) %>%
select(Fold, Repeat, rowIndex) %>%
inner_join(pred_ensemble_mediana, by = c("Repeat", "rowIndex"), suffix = c("","_") ) %>%
group_by(Fold, Repeat) %>%
nest() %>%
mutate( roc = map(data, ret_roc), auc = unlist(map(data, ret_roc_auc)) )
resultado_ensemble_pra_mostrar <- tibble( "Média AUC" = mean(pred_ensemble_mediana_foldado$auc), "Desvio-padrão AUC" = sd(pred_ensemble_mediana_foldado$auc) )
resultado_ensemble_pra_mostrar## # A tibble: 1 x 2
## `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl>
## 1 0.992 0.00748
É possível propor uma forma de avaliação do impacto das variáveis indireta e agnóstica ao modelo, ou seja, que pode ser usada da mesma forma para qualquer modelo e .
Para avaliar o impacto de cada variável nos resultados de cada um dos 5 modelos, podemos gerar casos sintéticos e avaliar a probabilidade de esses casos serem malignos segundo cada modelo.
Dois grupos de casos sintéticos foram gerados.
No primeiro grupo, para cada variável explicativa \(x_i\), foram criados casos onde o valor de \(x_i\) varia em relação a seus quantis, enquanto as outras permanecem parada no valor das suas medianas.
No segundo grupo, para avaliar uma região onde os casos têm mais chance de ser malignos, as variáveis que não estão sendo avaliadas a cada gráfico são mantidas no quantil 0,65.
medias_treino <- treino %>%
select(-Diagnosis) %>%
summarise_all(mean) %>%
gather(variavel, media)
dps_treino <- treino %>%
select(-Diagnosis) %>%
summarise_all(sd) %>%
gather(variavel, dp)
mediana_treino <- treino %>%
select(-Diagnosis) %>%
summarise_all(median) %>%
gather(variavel, mediana)
variaveis <- medias_treino %>%
select(variavel)
caso_media <- medias_treino %>%
spread(variavel, media)
quantis <- treino %>%
select(-Diagnosis) %>%
gather(variavel, valor) %>%
crossing( tibble( q = seq(0,1, 0.05)) ) %>%
group_by(variavel,q ) %>%
summarise( valor = quantile(valor, mean(q))) ## `summarise()` regrouping output by 'variavel' (override with `.groups` argument)
casos <- variaveis %>%
crossing(variaveis %>% rename(fixa = variavel)) %>%
crossing(tibble( q = seq(0,1, 0.05)) ) %>%
arrange(q, variavel, fixa) %>%
mutate(caso = (row_number()-1) %/% nrow(variaveis) + 1) %>%
mutate(q_cada_variavel = if_else(variavel == fixa, q, 0.5 )) %>%
inner_join(quantis, by=c("q_cada_variavel" = "q", "fixa" = "variavel") ) %>%
select(-q_cada_variavel) %>%
spread(fixa, valor)
prev_log <- predict(model_logistic, casos, type = "prob") %>%
mutate(modelo = "Logistica") %>%
bind_cols(casos)
prev_ranger <- predict(model_ranger, casos, type = "prob") %>%
mutate(modelo = "Floresta") %>%
bind_cols(casos)
prev_pca <- predict(model_net_pca, casos, type = "prob") %>%
mutate(modelo = "Lasso-Ridge PCA") %>%
bind_cols(casos)
prev <-
bind_rows(prev_log, prev_ranger, prev_pca)
prev %>%
ggplot() +
geom_line( aes(x = q, y = M, color = modelo) ) +
facet_wrap( ~variavel) +
theme_minimal() +
theme(
aspect.ratio = .4,
strip.text = element_text(size = 5),
strip.background.y = element_rect(size = c(0.1,0.1)),
strip.switch.pad.grid = unit(.01,"cm"),
strip.switch.pad.wrap = unit(.01,"cm"),
axis.text = element_text(size = 4)
) +
theme(legend.position = "top" )casos <- variaveis %>%
crossing(variaveis %>% rename(fixa = variavel)) %>%
crossing(tibble( q = seq(0,1, 0.05)) ) %>%
arrange(q, variavel, fixa) %>%
mutate(caso = (row_number()-1) %/% nrow(variaveis) + 1) %>%
mutate(q_cada_variavel = if_else(variavel == fixa, q, 0.65 )) %>%
inner_join(quantis, by=c("q_cada_variavel" = "q", "fixa" = "variavel") ) %>%
select(-q_cada_variavel) %>%
spread(fixa, valor)
prev_log <- predict(model_logistic, casos, type = "prob") %>%
mutate(modelo = "Logistica") %>%
bind_cols(casos)
prev_ranger <- predict(model_ranger, casos, type = "prob") %>%
mutate(modelo = "Floresta") %>%
bind_cols(casos)
prev_pca <- predict(model_net_pca, casos, type = "prob") %>%
mutate(modelo = "Lasso-Ridge PCA") %>%
bind_cols(casos)
prev <-
bind_rows(prev_log, prev_ranger, prev_pca)
prev %>%
ggplot() +
geom_line( aes(x = q, y = M, color = modelo) ) +
facet_wrap( ~variavel) +
theme_minimal() +
theme(
aspect.ratio = .4,
strip.text = element_text(size = 5),
strip.background.y = element_rect(size = c(0.1,0.1)),
strip.switch.pad.grid = unit(.01,"cm"),
strip.switch.pad.wrap = unit(.01,"cm"),
axis.text = element_text(size = 4)
) +
theme(legend.position = "top" )roc_best_modelo <- function(modelo)
{
dados_roc <- modelo$pred %>%
semi_join(modelo$bestTune) %>%
mutate(obs = if_else(obs == "B", 0, 1))
dados_curva <- roc(dados_roc$obs, dados_roc$M)
tibble( sensitivities = dados_curva$sensitivities,
specificities= dados_curva$specificities,
thresholds = dados_curva$thresholds,
auc = dados_curva$auc
)
}
pred_best_modelo <- function(modelo)
{
modelo$pred %>%
semi_join(modelo$bestTune)
}
plota_roc_modelo <- function(modelo, prob, nome)
{
roc_modelo <- roc_best_modelo(modelo)
pred_modelo <- pred_best_modelo(modelo)
threshold_escolhido <- roc_modelo %>%
filter(thresholds >= prob) %>%
top_n(n = 1, wt = desc(thresholds ))
ggplot(pred_modelo ) +
geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs ) ) +
geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
}
plota_roc_ensemble <- function(modelo, prob, nome)
{
prob <- 0.3
nome <- "Ensemble"
dados_curva <- roc(pred_ensemble_mediana$obs, pred_ensemble_mediana$M)
roc_modelo <- tibble( sensitivities = dados_curva$sensitivities,
specificities= dados_curva$specificities,
thresholds = dados_curva$thresholds,
auc = dados_curva$auc
)
pred_modelo <- pred_ensemble_mediana
threshold_escolhido <- roc_modelo %>%
filter(thresholds >= prob) %>%
top_n(n = 1, wt = desc(thresholds ))
ggplot(pred_modelo ) +
geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs ) ) +
geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
}
plota_roc_modelo(model_logistic, 1e-8, "Logistica")plota_roc_ensemble <- function(modelo, prob, nome)
{
prob <- 0.17
nome <- "Ensemble"
dados_curva <- roc(pred_ensemble_mediana$obs, pred_ensemble_mediana$M)
roc_modelo <- tibble( sensitivities = dados_curva$sensitivities,
specificities= dados_curva$specificities,
thresholds = dados_curva$thresholds,
auc = dados_curva$auc
)
pred_modelo <- pred_ensemble_mediana
threshold_escolhido <- roc_modelo %>%
filter(thresholds >= prob) %>%
top_n(n = 1, wt = desc(thresholds ))
ggplot(pred_modelo ) +
geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs ) ) +
geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
}
plota_roc_ensemble()confusao_teste <- function(modelo, threshold)
{
prev_teste_svm <- predict(modelo, teste, type = "prob")
pred_teste_svm <- prev_teste_svm %>%
mutate(pred = as.factor( if_else(M > threshold, "M", "B"))) %>%
select(pred)
conf_svm <- confusionMatrix(data = pred_teste_svm$pred, reference = teste$Diagnosis, positive = "M")
}
gera_dados_confusao_teste <- function(modelo, threshold)
{
prev_teste_svm <- predict(modelo, teste, type = "prob")
pred_teste_svm <- prev_teste_svm %>%
mutate(pred = as.factor( if_else(M > threshold, "M", "B"))) %>%
select(pred)
#tibble(data = pred_teste_svm$pred, reference = teste$Diagnosis, positive = "M")
bind_cols(pred_teste_svm, teste)
}## Confusion Matrix and Statistics
##
## Reference
## Prediction B M
## B 89 7
## M 10 65
##
## Accuracy : 0,9006
## 95% CI : (0,8456, 0,941)
## No Information Rate : 0,5789
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0,7972
##
## Mcnemar's Test P-Value : 0,6276
##
## Sensitivity : 0,9028
## Specificity : 0,8990
## Pos Pred Value : 0,8667
## Neg Pred Value : 0,9271
## Prevalence : 0,4211
## Detection Rate : 0,3801
## Detection Prevalence : 0,4386
## Balanced Accuracy : 0,9009
##
## 'Positive' Class : M
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction B M
## B 91 1
## M 8 71
##
## Accuracy : 0,9474
## 95% CI : (0,9024, 0,9757)
## No Information Rate : 0,5789
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0,8935
##
## Mcnemar's Test P-Value : 0,0455
##
## Sensitivity : 0,9861
## Specificity : 0,9192
## Pos Pred Value : 0,8987
## Neg Pred Value : 0,9891
## Prevalence : 0,4211
## Detection Rate : 0,4152
## Detection Prevalence : 0,4620
## Balanced Accuracy : 0,9527
##
## 'Positive' Class : M
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction B M
## B 91 2
## M 8 70
##
## Accuracy : 0,9415
## 95% CI : (0,8951, 0,9716)
## No Information Rate : 0,5789
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0,8814
##
## Mcnemar's Test P-Value : 0,1138
##
## Sensitivity : 0,9722
## Specificity : 0,9192
## Pos Pred Value : 0,8974
## Neg Pred Value : 0,9785
## Prevalence : 0,4211
## Detection Rate : 0,4094
## Detection Prevalence : 0,4561
## Balanced Accuracy : 0,9457
##
## 'Positive' Class : M
##
pred_test_1 <- predict(model_logistic, teste, type = "prob") %>%
mutate(row = row_number())
pred_test_2 <- predict(model_net_pca, teste, type = "prob") %>%
mutate(row = row_number())
pred_test_3 <- predict(model_ranger, teste, type = "prob") %>%
mutate(row = row_number())
pred_teste_ensemble <- bind_rows(pred_test_1, pred_test_2, pred_test_3) %>%
group_by(row) %>%
summarise (M = median(M)) %>%
mutate(pred = as.factor( if_else(M > 0.17, "M", "B"))) ## `summarise()` ungrouping output (override with `.groups` argument)
confusao_ensemble <- confusionMatrix(data = pred_teste_ensemble$pred, reference = teste$Diagnosis, positive = "M")
confusao_ensemble ## Confusion Matrix and Statistics
##
## Reference
## Prediction B M
## B 91 1
## M 8 71
##
## Accuracy : 0,9474
## 95% CI : (0,9024, 0,9757)
## No Information Rate : 0,5789
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0,8935
##
## Mcnemar's Test P-Value : 0,0455
##
## Sensitivity : 0,9861
## Specificity : 0,9192
## Pos Pred Value : 0,8987
## Neg Pred Value : 0,9891
## Prevalence : 0,4211
## Detection Rate : 0,4152
## Detection Prevalence : 0,4620
## Balanced Accuracy : 0,9527
##
## 'Positive' Class : M
##
dados_confusao_teste1 <- gera_dados_confusao_teste(model_logistic, 1e-8) %>%
mutate(modelo = "Logistic", linha = row_number())
dados_confusao_teste2 <- gera_dados_confusao_teste(model_net_pca, 0.18) %>%
mutate(modelo = "Lasso-Ridge", linha = row_number())
dados_confusao_teste3 <- gera_dados_confusao_teste(model_ranger , 0.25) %>%
mutate(modelo = "Random Forest", linha = row_number())
dados_confusao_teste_todos <- bind_rows(
dados_confusao_teste1,
dados_confusao_teste2,
dados_confusao_teste3
)
pred_errados_teste <- dados_confusao_teste_todos %>%
filter(pred != Diagnosis) %>%
group_by(linha, pred, Diagnosis) %>%
summarise(n= n()) %>%
arrange(desc(n))## `summarise()` regrouping output by 'linha', 'pred' (override with `.groups` argument)
## # A tibble: 23 x 4
## # Groups: linha, pred [23]
## linha pred Diagnosis n
## <int> <fct> <fct> <int>
## 1 37 M B 3
## 2 47 M B 3
## 3 111 B M 3
## 4 148 M B 3
## 5 150 M B 3
## 6 94 M B 2
## 7 120 M B 2
## 8 139 M B 2
## 9 2 B M 1
## 10 10 M B 1
## # ... with 13 more rows
dados_confusao_teste_todos_tidy <- dados_confusao_teste_todos %>%
gather(atributo, valor, - pred,-Diagnosis, - modelo, - linha)
dados_confusao_teste_todos_certos <-
dados_confusao_teste_todos_tidy %>%
filter(pred == Diagnosis)
dados_confusao_teste_todos_falso_negativo <-
dados_confusao_teste_todos_tidy %>%
filter(pred == "M" & Diagnosis == "B")
dados_confusao_teste_todos_falso_positivo <-
dados_confusao_teste_todos_tidy %>%
filter(pred == "B" & Diagnosis == "M")
ggplot() +
geom_jitter(data = dados_confusao_teste_todos_certos %>% filter(Diagnosis == "M"), aes(x = valor, y = 0), alpha = 0.05, size = 0.5, color = "red") +
geom_jitter(data = dados_confusao_teste_todos_certos %>% filter(Diagnosis == "B") , aes(x = valor, y = 0), alpha = 0.05, size = 0.5, color = "blue") +
geom_jitter(data = dados_confusao_teste_todos_falso_negativo, aes(x = valor, y = 0), alpha = 0.5, color = "blue", size = 0.5) +
geom_jitter(data = dados_confusao_teste_todos_falso_positivo , aes(x = valor, y = 0), alpha = 0.5, color = "red", size = 0.5) +
facet_wrap(~atributo, scales = "free") +
theme_minimal() +
theme(
aspect.ratio = .4,
strip.text = element_text(size = 5),
strip.background.y = element_rect(size = c(0.1,0.1)),
strip.switch.pad.grid = unit(.01,"cm"),
strip.switch.pad.wrap = unit(.01,"cm"),
axis.text = element_text(size = 4) ,
axis.text.y = element_blank()
)Shiny é um pacote R que permite a criação de aplicações WEB interativas que ativam código em R, com as possibilidades de formatação de tabelas, geração de gráficos, execução de modelos etc que as bibliotecas do R disponibilizam.
Os códigos desta seção que fala do Shiny não são executados diretamente, apenas são mostrados.
Uma aplicação Shiny mínima tem 4 itens:
Carregamento da biblioteca shiny
Definição da interface do usuário, ou seja, a página HTML com a qual o usuário vai interagir. Essa interface é definida com o uso da função fluidPage()
Definição do comportamento da aplicação com a função server()
Execução da função shinyApp(ui, server) para construir e iniciar uma aplicação.
Abaixo uma aplicação Shiny simples
Na função ui() do aplicativo shiny, definimos a interface.
No exemplo abaixo, além dos títulos, há um select e espaços para um gráfico e duas tabelas
Na função server() do aplicativo shiny, definimos o comportamento
No exemplo abaixo, definimos como o gráfico e as tabelas são gerados a partir do conteúdo do select contido na interface.
Repare que há duas tabelas, cada uma gerada com uma biblioteca diferente.
A biblioteca reactable permite que tenhamos mais controle sobre a estética
A biblioteca rhandontable parece mais com uma tabela de Excel, permitindo inclusive copiar e colar os números.
selectPaises <-
selectInput(
"pais",
label = "País",
choices = gapminder$country,
multiple = TRUE
)
ui <- fluidPage(
h1("Dashboard Gapminder"),
hr(),
h3("Filtros:"),
selectPaises,
h3("Gráfico:"),
plotOutput("grafico"),
h3("Tabela bonita:"),
reactableOutput("tabela_reactable"),
h3("Tabela pros viciados:"),
rHandsontableOutput("tabela_rhanson")
)
server <- function(input, output, session) {
output$grafico <- renderPlot({
gapminder %>%
filter(country %in% input$pais) %>%
ggplot() +
geom_line(aes(x = year, y = gdpPercap, color = country )) +
geom_point(aes(x = year, y = gdpPercap, color = country )) +
theme_light()
})
output$tabela_reactable <- renderReactable({
gapminder %>%
filter(country %in% input$pais) %>%
select(
country,
year,
gdpPercap
) %>%
pivot_wider(
names_from = year,
values_from = gdpPercap
) %>%
reactable(
defaultColDef = colDef(
format = colFormat(digits = 0, separators = TRUE)
)
)
})
output$tabela_rhanson <- renderRHandsontable({
gapminder %>%
filter(country %in% input$pais) %>%
select(
country,
year,
gdpPercap
) %>%
pivot_wider(
names_from = year,
values_from = gdpPercap
) %>%
mutate_at(
vars(matches("[0-9]{4}")),
~round(x = .x, digits = 0)
) %>%
rhandsontable(
readOnly = TRUE
) %>%
hot_cols(
format = "0,000",
language = "pt-BR"
)
})
}
shinyApp(ui, server)Diferentemente da programação feita em scripts, a execução no Shiny não segue a ordem dos comandos. É o conceito da programação reativa
A execução das funções do tipo render (que já vimos), reactive, observe e observeEvent é preguiçosa. Só acontece quando os inputs são alterados.
No caso do nosso aplicativo simples, apenas quando o input “pais” é alterado, as funções render relativas aos outputs (gráfico e tabela) são executadas.
Temos 2 papéis para os objetos: produtores e consumidores
Os inputs são produtores e os outputs são consumidores.
Vamos começar a ver agora objetos que funcionam tanto como produtores como consumidores.
As reactive expressions do shiny são objetos que assumem os dois papéis.
Elas possibilitam que cálculos e tratamentos que são comuns a várias saídas sejam executados apenas uma vez.
Na aplicação que desenhamos anteriormente, o select de país filtra da mesma forma os dados que são mostrados no gráfico e nas tabelas.
Podemos inserir um reactive que vai realizar esse tratamento que é comum a todas as saídas de uma só vez.
Assim evitamos mais facilmente, também, duplicar código.
Esta é a sintaxe das reactive expressions.
selectPaises <-
selectInput(
"pais",
label = "País",
choices = gapminder$country,
multiple = TRUE
)
ui <- fluidPage(
h1("Dashboard Gapminder"),
hr(),
h3("Filtros:"),
selectPaises,
h3("Gráfico:"),
plotOutput("grafico"),
h3("Tabela bonita:"),
reactableOutput("tabela_reactable"),
h3("Tabela pros viciados:"),
rHandsontableOutput("tabela_rhanson")
)
server <- function(input, output, session) {
dados_filtrados <- reactive({
gapminder %>%
filter(country %in% input$pais)
})
output$grafico <- renderPlot({
dados_filtrados() %>%
ggplot() +
geom_line(aes(x = year, y = gdpPercap, color = country )) +
geom_point(aes(x = year, y = gdpPercap, color = country )) +
theme_light()
})
output$tabela_reactable <- renderReactable({
dados_filtrados() %>%
select(
country,
year,
gdpPercap
) %>%
pivot_wider(
names_from = year,
values_from = gdpPercap
) %>%
reactable(
defaultColDef = colDef(
format = colFormat(digits = 0, separators = TRUE)
)
)
})
output$tabela_rhanson <- renderRHandsontable({
dados_filtrados() %>%
select(
country,
year,
gdpPercap
) %>%
pivot_wider(
names_from = year,
values_from = gdpPercap
) %>%
mutate_at(
vars(matches("[0-9]{4}")),
~round(x = .x, digits = 0)
) %>%
rhandsontable(
readOnly = TRUE
) %>%
hot_cols(
format = "0,000",
language = "pt-BR"
)
})
}
shinyApp(ui, server)Este dashboard de exemplo tem algumas funcionalidades interessantes que vamos explorar nos próximos slides
library(shiny)
library(tidyverse)
library(reactable)
library(sf)
library(waiter)
library(ggmap)
library(RColorBrewer)
library(scales)
#### LEITURA DE DADOS ####
gd <- read_rds("dados/gd/gd.rds") %>%
mutate(
unidade = 1
)
#### VARIÁVEIS DE ESTADO GLOBAIS ####
ultimos_pontos <- gd
ultimo_grafico <- NA
ultimo_n_amostra <- 0
#### OPÇÕES DOS SELECTS ####
agrupadores <- c(
"Setor econômico" = "USER_DscCl",
"Tarifa" = "USER_DscGr",
"Tipo de geração" = "USER_DscMo",
"Município" = "USER_NomMu",
"Região" = "USER_NomRe",
"UF" = "USER_SigUF",
"Tipo de fonte" = "USER_SigTi",
"Combustível" = "USER_DscCo",
"Agente" = "USER_SigAg"
)
somas <- c(
"Quantidade" = "unidade",
"Potência (MW)" = "USER_MdaPo"
)
#### ELEMENTOS DE UI ####
mapa <- plotOutput(
"mapa",
height = "600px",
brush = brushOpts(
id = "brush_mapa",
#resetOnNew = TRUE,
opacity = 0.2,
stroke = "black",
fill = "white",
clip = FALSE,
delay = 5000
)
)
selectUFs <- selectInput(
"ufs",
"UFs:",
choices = unique(gd$USER_SigUF) %>% sort(),
multiple = TRUE
)
slider_n_amostra <- sliderInput(
"n_amostra",
label = "% amostrado no gráfico",
min = 10,
max = 100,
step = 10,
value = 10,
post = "%"
)
slider_tamanho_ponto <- sliderInput(
"tamanho_ponto",
"Tamanho do ponto",
min = 0,
max = 2,
step = 0.05,
value = 0.05
)
agrupador <- selectInput(
"campo_grupo",
"Agrupar por:",
choices = agrupadores
)
soma_por <- selectInput(
"campo_soma",
"Somar por:",
choices = somas
)
tabela_count <- reactableOutput("info")
grafico_barras <- plotOutput("grafico_barras")
#### DIAGRAMAÇÃO DA UI ####
ui <- fluidPage(
theme = "bootstrap.css",
use_waiter(),
titlePanel(
fluidRow(
column(width = 1, img(height = 40, src = "logo-epe-azul-15-anos.gif")),
column(
width = 11,
tags$div(style = "height:10px"),
h3("Dashboard Geração Distribuída")
)
)
),
tags$hr(),
sidebarLayout(
sidebarPanel(
h4("Filtros" %>% tags$b()),
wellPanel(
selectUFs
),
h4("Configurações do gráfico" %>% tags$b()),
wellPanel(
agrupador,
soma_por,
slider_n_amostra,
slider_tamanho_ponto
),
width = 3
),
mainPanel(
width = 9,
splitLayout(
cellWidths = c("60%","40%"),
div(style = 'overflow: hidden',mapa),
div(style = 'overflow: hidden',
verticalLayout(
grafico_barras,
tabela_count
)
)
)
)
)
)
server <- function(input, output, session) {
#### MAPA ####
output$mapa <- renderPlot({
pontos_selecionados <- pontos()
if(isTRUE(all_equal(ultimos_pontos, pontos_selecionados)) & ultimo_n_amostra == input$n_amostra)
{
resposta <- ultimo_grafico
}
else
{
print("atualiza")
xmin <- min(pontos_selecionados$X)
ymin <- min(pontos_selecionados$Y)
xmax <- max(pontos_selecionados$X)
ymax <- max(pontos_selecionados$Y)
margem_x <- min((ymax - ymin)/7,2)
margem_y <- min((xmax - xmin)/7,2)
xmin <- xmin - margem_x
ymin <- ymin - margem_y
xmax <- xmax + margem_x
ymax <- ymax + margem_y
if (is.na(xmin)){
local <- c(-67.85551, -31.76278, -34.80921, -1.198088)
}
else{
local <- c(xmin, ymin, xmax, ymax)
}
map = get_map(local,source = "stamen", maptype = "toner-lite")
resposta <- ggmap(ggmap = map) +
stat_density_2d(
aes(x = X, y = Y, fill = ..level.. ),
bins = 30,
geom = "polygon",
data = pontos_selecionados ,
alpha = .1
) +
geom_point(
data = pontos_selecionados,
aes(
x = X,
y = Y
),
alpha = 0.01,
color = "darkblue",
size = input$tamanho_ponto
) +
scale_fill_gradientn(colors = brewer.pal(7, "YlOrRd")) +
theme_inset()
ultimo_grafico <<- resposta
ultimos_pontos <<- pontos_selecionados
ultimo_n_amostra <<- input$n_amostra
}
resposta
})
#### PONTOS DO MAPA SELECIONADOS ####
pontos <- reactive({
waiter <- Waiter$new(id = "mapa")$show()
resposta <- brushedPoints(gd_sample(), input$brush_mapa, xvar = "X", yvar = "Y")
if (nrow(resposta) == 0){
resposta <- gd_sample()
}
resposta
})
#### DADOS SELECIONADOS PELO MAPA####
dados <- reactive({
resposta <- brushedPoints(gd_filtro(), input$brush_mapa, xvar = "X", yvar = "Y")
if (nrow(resposta) == 0){
resposta <- gd_filtro()
}
campo_grupo <- names(agrupadores)[agrupadores == input$campo_grupo]
campo_soma <- names(somas)[somas == input$campo_soma]
resposta %>%
group_by_at(input$campo_grupo) %>%
summarise_at(input$campo_soma, sum) %>%
ungroup() %>%
mutate(
Parcela = .data[[input$campo_soma]]/sum(.data[[input$campo_soma]])
) %>%
arrange(desc(.data[[input$campo_soma]])) %>%
rename(
!!campo_grupo := 1,
!!campo_soma := 2
)
})
#### TABELA ####
output$info <- renderReactable({
dados() %>%
reactable(
pagination = FALSE,
defaultColDef = colDef(
format = colFormat(
digits = 0,
separators = TRUE
)
),
columns = list(
Parcela = colDef(
format = colFormat(
percent = TRUE,
digits = 0
)
)
)
)
})
#### GRÁFICO DE BARRAS HORIZONTAIS ####
output$grafico_barras <- renderPlot({
campo_grupo <- names(agrupadores)[agrupadores == input$campo_grupo]
campo_soma <- names(somas)[somas == input$campo_soma]
dados() %>%
mutate(
!!campo_grupo :=
fct_lump(
f = .data[[campo_grupo]],
prop = 0.05,
w = .data[[campo_soma]] ,
other_level = "Outros"
)
) %>%
mutate(
!!campo_grupo :=
fct_reorder(
.f = .data[[campo_grupo]],
.x = .data[[campo_soma]],
.fun = sum
)
) %>%
ggplot() +
geom_col(
aes(
x = .data[[campo_grupo]],
y = .data[[campo_soma]],
fill = .data[[campo_grupo]]
)
) +
labs(
y = campo_soma,
x = campo_grupo,
fill = campo_grupo
) +
coord_flip() +
theme_minimal() +
scale_y_continuous(
labels = number_format(
big.mark = ".",
accuracy = 1
)
) +
theme(
legend.position = "top"
) +
NULL
})
gd_filtro <- reactive({
gd %>%
filter(USER_SigUF %in% input$ufs | length(input$ufs) == 0 )
})
gd_sample <- reactive({
gd_filtro() %>%
sample_frac(input$n_amostra/100)
})
}
shinyApp(ui, server)O usuário pode selecionar áreas do mapa.
A seleção dessa área leva a uma seleção de pontos
Essa opção é habilitada com brushOpts e os pontos selecionados são acessados com brushedPoints
Este tipo de interação pode ser usado com qualquer tipo de gráfico
Esta aplicação lida com mais de 100.000 dados a respeito de instalações de geração distribuída.
Portanto, o gráfico demora um pouco para ser plotado.
Para que o usuário saiba que a aplicação está calculando alguma coisa e não travada, usamos um waiter
Algumas variáveis de estado são criadas para ajudar a controlar a aplicação.
Elas são atribuídas com o operador <<-
O Tidy evaluation deixa que usemos os nomes dos campos sem aspas.
Em aplicações Shiny exploratórias, é comum que os campos que participam dos agrupamentos ou dos cálculos sejam escolhidos pelo usuário. Portanto, estes campos estarão no comteúdo de variáveis.
Para que as funções do tidyverse usem o conteúdo da variável como o campo a ser considerado (e não o próprio nome da variável), temos que lançar mão das funções group_by_at, summarise_at, .data[[variável]] e o operador !!.
A diagramação da página é toda feita com as funções de diagramação do Shiny.
Entretanto, cada objeto é definido anteriormente, o que é diferente do template que é gerado pelo RStudio quando um template Shiny é gerado.
Isso parece besteira, mas separar a criação dos objetos da diagramação destes objetos facilita muito uma reorganizaçào da diagramação e causa menos problemas relacionados a erros de sintaxe como falta de vírgulas e parênteses
Os pacotes integrantes do conjunto tidyverts tratam as séries temporais do jeito tidy
| -Funções básicas para o tsibble (tibble temporal) |
-Funções pra extrações de features e estatísticas |
-Funções para projeção |
“)tsibbler a_r(” é uma classe derivada da tibbl que tem funcionalidades para séries temporais
A função “)as_tsibble()r a_r(” cria um tibble a partir de vários tipos de objetos que armazenam séries temporais
O tsibble representa bem séries temporais regulares e reconhece o intervalo entre os valores
## # A tsibble: 492 x 2 [1M]
## index value
## <mth> <dbl>
## 1 1980 jan 0.0662
## 2 1980 fev 0.0462
## 3 1980 mar 0.0604
## 4 1980 abr 0.0529
## 5 1980 mai 0.057
## 6 1980 jun 0.0531
## 7 1980 jul 0.0555
## 8 1980 ago 0.0495
## 9 1980 set 0.0423
## 10 1980 out 0.0948
## # ... with 482 more rows
No código abaixo calculamos o número de instalações de Geração Distribuída por UF as partir dos dados a respeito das instalações.
Repare que a função ")`yearmonth()### cria uma data mensal que é específica para períodos de um mês
dados_gd <- read_rds("dados/gd/gd.rds")
meses <- tibble(mes = seq.Date(from = make_date(2009,1,1), to = make_date(2020,3,1), by = "month" ) %>% yearmonth())
dados_gd_cum <- dados_gd %>%
rename(
uf = USER_SigUF
) %>%
crossing(meses) %>%
filter(DthConexao < mes) %>%
group_by(
uf,
mes
) %>%
summarise(
n = n(),
potencia = sum(USER_MdaPo)
)## Warning in mask$eval_all_filter(dots, env_filter): Métodos incompatíveis
## ("<.Date", "<.vctrs_vctr") para "<"
## `summarise()` regrouping output by 'uf' (override with `.groups` argument)
A função as_tibble aceita um campo como referência da data, via parâmetro index, e um campo como chave de cada série contida no tibble: key
## # A tsibble: 2.179 x 4 [1M]
## # Key: uf [28]
## # Groups: uf [28]
## uf mes n potencia
## <chr> <mth> <int> <dbl>
## 1 AC 2015 mai 1 4
## 2 AC 2015 jun 1 4
## 3 AC 2015 jul 2 6
## 4 AC 2015 ago 3 38.5
## 5 AC 2015 set 3 38.5
## 6 AC 2015 out 3 38.5
## 7 AC 2015 nov 3 38.5
## 8 AC 2015 dez 3 38.5
## 9 AC 2016 jan 3 38.5
## 10 AC 2016 fev 3 38.5
## # ... with 2.169 more rows
Podemos mudar a periodicidade da série indexando por uma transformação do índice atual, usando ")`group_by_key** e index_by.
No caso, queremos as informações no final do ano, portanto vamos usar **last()`r a_r(".
gd_anual <- dados_gd_tsbl %>%
group_by_key() %>%
index_by(ano = year(mes) ) %>%
summarise(
n = last(n),
potencia = last(potencia)
)
gd_trimestre <- dados_gd_tsbl %>%
group_by_key() %>%
index_by(trimestre = yearquarter(mes)) %>%
summarise(
n = last(n),
potencia = last(potencia)
)
gd_anual## # A tsibble: 214 x 4 [1Y]
## # Key: uf [28]
## uf ano n potencia
## <chr> <dbl> <int> <dbl>
## 1 AC 2015 3 38.5
## 2 AC 2016 4 40.5
## 3 AC 2017 12 173.
## 4 AC 2018 38 630.
## 5 AC 2019 134 1968.
## 6 AC 2020 172 2353.
## 7 AL 2015 3 57.1
## 8 AL 2016 18 270.
## 9 AL 2017 77 1195.
## 10 AL 2018 255 3203.
## # ... with 204 more rows
## # A tsibble: 737 x 4 [1Q]
## # Key: uf [28]
## uf trimestre n potencia
## <chr> <qtr> <int> <dbl>
## 1 AC 2015 Q2 1 4
## 2 AC 2015 Q3 3 38.5
## 3 AC 2015 Q4 3 38.5
## 4 AC 2016 Q1 3 38.5
## 5 AC 2016 Q2 3 38.5
## 6 AC 2016 Q3 4 40.5
## 7 AC 2016 Q4 4 40.5
## 8 AC 2017 Q1 4 40.5
## 9 AC 2017 Q2 7 75.8
## 10 AC 2017 Q3 9 121.
## # ... with 727 more rows
A função ")`autoplot()### cuida das principais tarefas para plotar uma série
#
#
# tsbl_ipca %>%
# filter(index > make_date(1994,7,1)) %>%
# autoplot(value, color = "darkblue", size = 1) +
# scale_y_continuous(labels = percent_format(accuracy = 1) ) +
# scale_x_date(date_breaks = "1 year", labels = year) +
# labs(
# x = "",
# y = "IPCA"
# ) +
# theme_minimal() +
# theme(
# axis.text.x = element_text(angle = 90)
# )
# A função ")gg_season") ajuda a fazer um gráfico mostrando os períodos sazonais para comparação
# tsbl_ipca %>%
# filter(year(index) > 2009) %>%
# gg_season(value, period = "year", labels = "both") +
# scale_y_continuous(
# labels = percent_format(accuracy = 0.1),
# breaks = seq(0.000, by = 0.001,to = 0.015) ) +
# scale_x_date(
# date_labels = "%b",
# date_breaks = "1 month"
# ) +
# labs(
# x = "",
# y = "IPCA"
# ) +
# theme_minimal() +
# theme(
# axis.text.x = element_text(angle = 90)
# )
#
#
#
# # tsbl_ipca %>%
# filter(year(index) > 1995) %>%
# gg_season(value, period = "year") +
# scale_y_continuous(
# labels = percent_format(accuracy = 0.1),
# breaks = seq(0.000, by = 0.002,to = 0.03) ) +
# scale_x_date(
# date_labels = "%b",
# date_breaks = "1 month"
# ) +
# labs(
# x = "",
# y = "IPCA"
# ) +
# theme_minimal() +
# theme(
# axis.text.x = element_text(angle = 90)
# )
# Tudo que pode ser automatizado deve ser automatizado devtools
Pacotes a instalar para começar:
Para checar se o ambiente está preparado para criar pacotes:
New Project -> New Directory -> R Package
The main difference between library() and require() is what happens if a package isn’t found. While library() throws an error, require() prints a warning message and returns FALSE. In practice, this distinction isn’t important because when building a package you should NEVER use either inside a package. See package dependencies for what you should do instead.
[.
No caso, queremos as informações no final do ano, portanto vamos usar ](#(257))
antipadrões de visualização e dados
dicas para boa visualização de dados
execução em paralelo no exemplo
função de probablidade acumulada empírica
lead() no contexto da group_by
PCA (Principal Component Analysis)
position - parâmetro da geom_col()
programação funcional no exemplo
seleção negativa de elementos com -
slice_max() no contexto da group_by()